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Abstract 


Objectives 

We  aimed  to  develop  and  evaluate  practical  methods  to  estimate  species  richness  and 
occupancy  of  diverse  taxonomic  groups  across  space  and  time  and  in  diverse  ecosystems.  As 
reflected  in  the  statement  of  need  to  which  this  project  responded,  our  work  is  relevant  to  the 
needs  of  the  Department  of  Defense  (DoD)  to  assess  and  monitor  native  species  and  to  evaluate 
the  potential  effects  of  land  use  and  management  actions  on  native  species.  Species  richness,  or 
the  number  of  native  species,  is  a  common  surrogate  measure  of  ecosystem  status.  Occupancy, 
the  probability  that  a  given  location  is  occupied  by  a  species,  can  serve  as  a  surrogate  measure  of 
the  species’  abundance,  which  in  turn  is  related  to  probability  of  persistence. 

This  project  addressed  five  major  objectives  consistent  with  the  primary  aim.  First,  we 
assessed  relations  between  environmental  variables  and  species  richness  at  nested  spatial  extents. 
Second,  we  developed  guidelines  for  periodicity  of  sampling.  Third,  we  tested  methods  for 
estimating  species  richness  of  multiple  taxonomic  groups  on  the  basis  of  spatial  variation  in 
occurrence  or  heterogeneity  of  one  group.  Fourth,  we  examined  the  extent  to  which  occupancy 
could  provide  a  foundation  for  measuring  species  richness.  Fifth,  we  investigated  potential 
responses  of  species  richness  and  occupancy  to  manageable  environmental  changes,  and 
meaningful  scales  of  sampling  for  detection  of  biological  effects  of  environmental  change. 

Technical  approach 

As  land  use  changes  in  the  multi -jurisdictional  lands  surrounding  DoD  installations  in  the 
United  States,  it  becomes  more  challenging  to  sustain  military  readiness  and  conserve  biological 
diversity.  We  sampled  anurans,  birds,  and  butterflies  in  the  Chesapeake  Bay  Lowlands 
(including  Joint  Base  Langley-Eustis),  and  birds  and  butterflies  in  the  central  and  western  Great 
Basin  (including  Hawthorne  Army  Depot  and  the  Marine  Corps  Mountain  Warfare  Training 
Center),  from  2012-2015.  Working  in  ecosystems  in  which  land-cover  configurations  and 
drivers  of  those  configurations  differ,  and  on  disparate  faunal  groups,  allowed  us  to  evaluate  the 
geographic  and  taxonomic  transferability  of  our  methods  and  inferences.  We  also  capitalized  on 
existing  data  and  local  experience  that  served  the  project  aims. 

Our  sampling  methods  allowed  us  to  estimate  detection  probability — the  probability  of 
detecting  a  species  if  it  is  present.  If  detection  probability  is  not  quantified,  inferred  relations 
between  environmental  covariates  and  species  occurrence  or  demography  may  be  erroneous.  We 
conducted  hierarchical  analyses  of  the  species  richness  of  birds  and  butterflies  in  the  central  and 
western  Great  Basin.  We  modeled  species  richness  at  two  nested  spatial  extents  as  functions  of 
environmental  variables  measured  at  either  extent.  We  used  the  same  data  to  develop  a  method  to 
identify  indicator  species  (small  sets  of  species  with  occurrence  patterns  that  are  related  to 
species  richness  of  larger  sets)  and  environmental  variables  that  explain  considerable  variation  in 
species  richness.  We  conducted  novel,  rigorous  external  evaluations  of  the  spatial  and  temporal 
transferability  of  such  indicator  species.  Furthermore,  we  evaluated  the  extent  to  which  different 
ecological  processes  might  drive  spatial  and  temporal  variation  in  species  identities. 

We  examined  the  responses  of  estimated  detection  probability  and  occupancy  of  birds  in 
the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin  to  the  duration  of  sampling. 
Additionally,  we  examined  whether  detection  probability  and  occupancy  of  birds  and  anurans  in 
the  Chesapeake  Bay  Lowlands,  and  butterflies  in  the  Great  Basin,  differed  if  sites  were  sampled 
repeatedly  on  a  single  day  versus  once  each  on  multiple  days. 


We  explored  the  extent  to  which  occupancy  of  butterflies  in  the  Chesapeake  Bay 
Lowlands  and  central  and  western  Great  Basin  was  associated  with  vegetation  structure  and 
composition,  topography,  and  other  environmental  attributes  and  whether  assumptions  of  closure 
were  met  with  assemblage-level  sampling  designs.  Furthermore,  we  identified  attributes  of 
wetlands  and  the  surrounding  terrestrial  environment  that  were  associated  with  anuran 
populations  of  high  relative  abundance  versus  populations  of  any  abundance  class,  and  evaluated 
how  these  attributes  relate  to  different  life  history  stages. 

We  examined  how  songbirds  in  the  Chesapeake  Bay  Lowlands,  including  species  of 
concern,  may  respond  to  changes  in  vegetation  fragmentation  and  structure  that  result  from 
urbanization  and  different  levels  of  use  by  white-tailed  deer.  We  also  explored  turnover  in 
species  composition  of  birds  and  butterflies  in  the  central  and  western  Great  Basin  to  gain  insight 
into  meaningful  scales  of  sampling  for  detection  of  biological  effects  of  environmental  change. 

Results 

In  the  central  and  western  Great  Basin,  species  richness  of  birds  and  butterflies  responded 
to  environmental  variables  at  different  spatial  extents.  Species’  identities  were  more  consistent  in 
space  in  butterfly  than  in  bird  assemblages,  whereas  spatial  nestedness  (the  extent  to  which 
species-poor  faunas  are  statistical  subsets  of  species-rich  faunas)  often  was  higher  in  butterfly 
than  in  bird  assemblages.  In  some  cases,  we  were  able  to  explain  more  than  80%  and  70%  of  the 
variation  in  species  richness  of  birds  and  butterflies,  respectively,  as  functions  of  vegetation  and 
topography.  The  predictive  accuracy  and  spatial  and  temporal  transferability  of  models  of 
species  richness  of  birds  and  butterflies  that  are  based  on  indicator  species  can  exceed  that  of 
models  based  on  environmental  variables. 

In  all  of  our  ecosystems,  modest  increases  in  the  duration  of  point  counts  for  birds  did  not 
affect  inferences  based  on  occupancy  models,  but  might  affect  estimates  of  species  richness  that 
were  not  detection-weighted.  Given  our  results,  we  believe  that  estimation  of  use  may  be  equally 
or  more  informative  than  estimation  of  occupancy  for  exploring  species-environment  relations. 

In  many  cases,  because  detection  probabilities  for  many  species  of  birds  and  butterflies 
were  low,  occupancy  could  not  be  estimated  for  more  than  50%  of  the  species  detected  in  each 
ecosystem  and  year.  Although  imperfect  detection  may  lead  to  biased  estimates  of  species 
occurrence,  demographic  parameters,  and  species-environment  relations,  achieving  certain 
scientific  objectives  still  may  require  use  of  occurrence  data  that  are  not  detection  weighted.  Our 
results  suggest  that  if  anurans  are  extirpated  from  wetlands  in  the  Chesapeake  Bay  Lowlands,  the 
wetlands  are  unlikely  to  be  recolonized. 

Benefits 

We  regularly  shared  data  and  inferences  with  environmental  managers  at  the  installations 
on  which  we  worked,  who  indicated  that  our  results  would  contribute  to  watershed  management 
programs  and  Integrated  Natural  Resources  Management  Plans.  We  also  shared  data  and 
inferences  with  DoD  partners,  including  Landscape  Conservation  Cooperatives.  In  the 
Chesapeake  Bay  Lowlands,  our  results  are  relevant  to  management  of  white-tailed  deer,  which 
have  substantial  effects  on  other  native  species  and,  via  transmission  of  disease  vectors,  public 
health.  In  the  Great  Basin,  our  results  are  applicable  to  vegetation  treatments  that  are  intended  to 
increase  the  probability  of  persistence  of  Greater  Sage-Grouse  and  benefit  hundreds  of  other 
species.  Our  results  also  suggest  directions  for  future  research  on  the  extent  to  which  simple 
measures  of  species  occurrence  may  provide  information  on  probabilities  of  species  persistence. 
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Objective 


The  desired  outcomes  of  the  statement  of  need  to  which  this  project  responded  were  knowledge 
that  improves  understanding  of  assessment,  monitoring,  and  management  of  native  species  at 
multiple  spatial  and  temporal  scales;  knowledge  that  can  assist  resource  managers  in  providing 
effective  and  cost-efficient  methods  for  accomplishing  such  assessment  and  monitoring  and  for 
projecting  the  effects  of  land  use  and  management  actions;  knowledge  that  addresses  diverse 
aspects  of  theory  related  to  assessment  and  monitoring  of  native  species  and  its  application 
(including  the  fact  that  habitat  is  species-specific,  and  therefore  different  species  will  have 
different  responses  to  environmental  change);  and  knowledge  that  promotes  assessment  and 
monitoring  of  native  species  within  and  among  land  management  and  administrative  units. 
Assessment  and  monitoring  of  native  species  is  relevant  to  the  Department  of  Defense  (DoD) 
because  maintenance  of  these  species  can  conflict  with  operations  that  sustain  military  readiness; 
because  the  DoD  is  committed  to  environmental  stewardship;  and  because  the  DoD  must  comply 
with  laws  that  protect  species  and  other  levels  of  biological  diversity.  Furthermore,  information 
on  where  and  when  species  of  concern  occur,  and  associations  between  those  species  and  then- 
environment,  often  is  necessary  to  inform  environmental  management  on  installations. 

In  response  to  the  statement  of  need,  the  objective  of  this  project  was  to  develop  and  validate 
practical  methods  to  estimate  species  richness  and  occupancy  across  space  and  time.  We  define 
species  richness  as  the  number  of  native  species  in  a  given  location  during  a  given  time  period. 
Occupancy,  the  probability  that  a  given  location  is  occupied  by  a  species,  is  a  state  variable  that 
can  serve  as  a  surrogate  measure  of  the  abundance  of  that  species  (MacKenzie  and  Nichols 
2004).  Abundance  and  geographic  distribution  are  strongly  related  to  probability  of  persistence 
(Lande  1993,  Foley  1994,  Harris  and  Pimm  2008).  We  aimed  to  develop  or  refine  methods  that 
facilitate  estimation  of  a  species’  spatial  pattern  of  occupancy,  detect  changes  in  its  distribution, 
and  make  inference  to  its  probability  of  regional  persistence.  We  also  aimed  to  develop  methods 
that  were  statistically  rigorous;  transferable  in  space  (including  within  and  among  ecosystems) 
and  time;  where  possible,  transferable  among  taxonomic  groups;  and  that  reconciled  statistical 
rigor  with  biological  realism.  We  quantified  the  extent  to  which  occupancy  of  three  faunal 
groups  in  three  ecosystems  (including  two  biogeographic  subregions  of  one  extensive  ecoregion) 
can  be  explained  as  a  function  of  geophysical  attributes,  land  cover,  and  the  composition  and 
structure  of  vegetation  at  multiple  spatial  extents  and  resolutions.  We  sampled  birds,  butterflies, 
and  anurans  in  the  Chesapeake  Bay  Lowlands  (including  Joint  Base  Langley-Eustis),  and  birds 
and  butterflies  in  the  central  and  western  Great  Basin  (including  Hawthorne  Army  Depot  and  the 
Marine  Corps  Mountain  Warfare  Training  Center).  We  are  continuing  to  use  the  results  to 
project  how  occupancy  and  species  richness  may  change  in  response  to  environmental  changes 
that  are  responsive  to  management.  We  are  providing  all  data  and  customized  summaries  to  DoD 
resource  managers  and  environmental  specialists,  and  installations  indicated  that  they  will 
incorporate  our  results  into  watershed  management  programs  and  Integrated  Natural  Resources 
Management  Plans. 

At  the  start  of  the  project,  we  had  five  major  tasks  that  responded  to  the  statement  of  need.  Our 
first  objective  was  to  assess  two  methods  for  correcting  bias  in  estimates  of  species  richness  at 
small  spatial  extents:  standardizing  sampling  effort  by  area  sampled  and  by  number  of 
individuals  sampled.  Over  the  course  of  the  project,  as  described  below,  we  refocused  this 
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objective  and  conducted  a  hierarchical  estimation  of  species  richness  of  birds  and  butterflies  in 
the  central  and  western  Great  Basin.  Our  second  objective  was  to  develop  guidelines  for 
periodicity  of  sampling.  We  used  data  on  anurans  from  the  Chesapeake  Bay  Lowlands  and  birds 
and  butterflies  from  the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin  to 
address  this  objective.  Third,  we  used  data  on  birds  and  butterflies  from  the  central  and  western 
Great  Basin  to  test  methods  for  estimating  species  richness  of  multiple  taxonomic  groups  on  the 
basis  of  spatial  variation  in  occurrence  or  heterogeneity  of  a  single  taxonomic  group.  Our  fourth 
objective  was  to  develop  a  structured  method  for  species-level  monitoring  that  used  occupancy 
as  a  surrogate  measure  of  abundance  for  estimating  patterns  of  species  richness.  We  used  data  on 
anurans  and  on  birds  and  butterflies  from  the  Chesapeake  Bay  Lowlands  and  central  and  western 
Great  Basin  to  address  that  objective.  Fifth,  we  aimed  to  model  potential  responses  of  species 
richness  and  occupancy  to  environmental  change.  In  the  Chesapeake  Bay  Lowlands,  we  focused 
on  responses  of  birds  to  urbanization  and  to  the  density  of  white-tailed  deer  ( Oclocoileus 
virginicinus ),  the  abundance  of  which  has  increased  by  more  than  3000%  over  the  past  century. 
Rapid  increases  in  the  regional  abundance  of  white-tailed  deer  have  been  linked  not  only  with 
changes  in  native  flora  and  fauna  but  with  increases  in  the  incidence  of  tick-borne  diseases  that 
affect  humans,  such  as  Lyme  disease.  The  abundance  and  density  of  deer  can  be  managed  by 
humans  both  within  and  beyond  installation  boundaries.  In  the  central  and  western  Great  Basin, 
our  results  are  applicable  to  planned  management  treatments  on  DoD  and  other  public  lands, 
such  as  manipulations  of  vegetation  that  are  intended  to  increase  the  probability  of  persistence  of 
Greater  Sage-Grouse  ( Centrocercus  urophasianus )  and  benefit  more  than  350  other  species  of 
animals  (BLM  2015,  USFS  2015). 
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Background 


Species  richness  commonly  is  used  as  a  surrogate  measure  of  the  status  of  an  ecological  system, 
and  it  is  a  measure  readily  understood  by  decision-makers  and  the  public.  Therefore,  species 
richness  long  has  been  a  focus  of  assessment  and  monitoring.  However,  species  richness  alone  is 
insufficient  to  characterize  whether  locations  are  likely  to  sustain  species  and  ecological  function 
over  the  long  term.  Additionally,  most  measures  of  species  richness  are  not  accurate  statistical 
estimates  because  it  is  impossible  to  determine  their  precision  or  the  effects  of  possible  biases, 
such  as  underestimation  of  species  richness  as  a  result  of  not  detecting  rare  or  cryptic  species. 
Nevertheless,  assessments  of  numerous  individual  species,  especially  when  density  or  abundance 
is  the  state  variable,  generally  have  been  considered  too  complex  and  expensive  to  be  feasible. 

Until  the  early  2000s,  most  methods  for  estimating  species  richness,  abundance,  and 
demographics  of  plant  and  animal  populations  did  not  account  for  detection  probability — the 
probability  of  detecting  a  species  at  a  given  location  conditional  on  its  presence.  If  detection 
probability  is  not  quantified,  estimators  of  presence,  abundance,  and  other  measures  of 
population  dynamics  may  be  biased  and  may  lead  to  erroneous  inferences  about  relations 
between  response  variables  and  environmental  covariates  (Gu  and  Swihart  2004,  MacKenzie 
2005).  Many  sampling  and  analytical  methods  have  been  developed  to  address  the  reality  of 
imperfect  detection  (e.g.,  Buckland  et  al.  2001,  Williams  et  al.  2002).  For  example,  the  single¬ 
season  occupancy  model  (MacKenzie  et  al.  2002)  was  developed  to  account  for  the  possibility  of 
false  absences — species  that  are  present,  but  not  detected,  during  a  given  sampling  period — in 
estimating  species’  spatial  or  temporal  distributions  (MacKenzie  et  al.  2006).  There  is  a  solid 
theoretical  basis  for  using  occupancy  as  a  surrogate  measure  of  a  species’  abundance  or  state 
(e.g.,  reproductive  status)  (He  and  Gaston  2003,  Nichols  et  al.  2007,  Green  et  al.  2011). 
Abundance,  in  turn,  is  strongly  related  to  probability  of  persistence  (Gaston  et  al.  2000, 
Zuckerberg  et  al.  2009),  and  maximizing  the  probability  of  species  persistence  commonly  is  a 
high  priority  for  federal  and  state  managers  of  natural  resources. 

Use  of  occupancy  models  has  been  advocated  not  only  on  the  basis  of  statistical  rigor,  but  on  the 
basis  of  feasibility.  The  data  required  to  estimate  demographic  parameters  or  abundance 
precisely,  especially  for  relatively  rare  species,  can  be  expensive  and  difficult  to  obtain 
(MacKenzie  et  al.  2004,  Corbani  et  al.  2014).  Because  occupancy  models  do  not  rely  on  captures 
or  detections  of  multiple  individuals  in  a  population,  they  can  be  applied  to  a  greater  proportion 
of  the  species  in  a  community.  Researchers  have  used  occupancy  models  to  estimate  species 
richness  (Kery  et  al.  2009,  Dorazio  et  al.  2010,  Carillo-Rubio  et  al.  2014),  trends  in  distribution 
(van  Strien  et  al.  2013,  Weir  et  al.  2014),  and  extinction-colonization  dynamics  (Fernandez- 
Chacon  et  al.  2014,  Pearson  et  al.  2015)  of  multiple  taxonomic  groups.  To  model  occupancy,  one 
selects  a  sample  of  sites  (e.g.,  quadrats  or  habitat  patches)  and  surveys  each  site  multiple  times  in 
a  given  season  or  year  (MacKenzie  et  al.  2002).  Multiple  surveys  are  necessary  to  estimate 
detection  probability,  and  the  number  and  timing  of  these  surveys  is  a  major  consideration  in  the 
study  design. 

There  are  four  core  assumptions  of  occupancy  models.  First,  occupancy  at  a  site  does  not  change 
among  surveys  (i.e.,  the  closure  assumption).  Second,  either  there  is  no  heterogeneity  in 
probabilities  of  detection  and  occupancy  of  a  given  species,  or  heterogeneity  is  effectively 
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modeled.  Third,  detections  of  a  given  species  are  independent  among  surveys  and  sites.  Fourth, 
the  species  of  interest  are  not  falsely  detected  (MacKenzie  et  al.  2006).  The  closure  assumption 
arguably  is  the  most  challenging.  The  occupancy  status  of  mobile  organisms  on  sample  units 
rarely  is  constant  across  the  season,  and  changes  in  occupancy  may  not  be  random.  Therefore, 
detection  estimates  that  treat  multiple  surveys  as  replicate  samples  of  the  same  species  may  be 
negatively  biased,  and  the  resulting  estimates  of  occupancy  positively  biased.  One  possibility  is 
to  sample  each  site  repeatedly  on  each  sampling  date  (the  robust  design;  Kendall  et  al.  1997). 
Another  option  is  to  relax  the  closure  assumption  by  assuming  that  species  are  available  for 
sampling  at  different  times  (Kendall  et  al.  2013).  We  explored  both  of  these  options  in  our  work. 
We  examined  trade-offs  between  single-day  and  multiple-days  sampling  of  anurans  and  birds  in 
the  Chesapeake  Bay  Lowlands  and  butterflies  in  the  central  and  western  Great  Basin,  and  we 
tested  whether  different  species  of  butterflies  in  our  three  ecosystems  met  the  assumption  of 
closure. 

Working  in  ecosystems  in  which  land-cover  configurations  and  drivers  of  those  configurations 
differ  allowed  us  to  evaluate  the  transferability  of  our  methods  and  inferences.  For  example, 
patches  of  deciduous  forest  in  the  Chesapeake  Bay  Lowlands  are  relatively  dispersed,  whereas 
patches  of  pinyon  and  juniper  woodland  in  the  Great  Basin  are  relatively  clumped;  riparian  land- 
cover  generally  is  more  extensive  and  less  fragmented  in  the  mesic  Chesapeake  Bay  Lowlands 
than  in  the  arid  Great  Basin.  Outside  of  installations,  which  have  relatively  contiguous  land 
cover,  most  of  the  Chesapeake  Bay  Lowlands  is  dominated  by  heavily  fragmented  private 
holdings,  whereas  the  Great  Basin  is  dominated  by  multijurisdictional  federal  lands. 

Throughout  the  mid-Atlantic  region,  including  the  Chesapeake  Bay  Lowlands,  fragmentation  of 
deciduous  forest  is  driven  by  urban,  exurban,  and  rural  settlement,  including  construction  of 
roads  and  other  infrastructure;  by  agriculture;  and  by  development  of  coal  and  wind  energy 
(Woolmer  et  al.  2008).  In  contrast,  fragmentation  of  woodlands  in  the  Intermountain  West, 
including  the  Great  Basin,  is  driven  both  by  natural  processes  such  as  wildfire  and  by  land  uses 
such  as  domestic  livestock  grazing,  energy  development  (wind,  solar,  oil,  gas,  and  thermal),  and 
exurbanization,  with  associated  infrastructure  such  as  secondary  roads.  Secondary  roads  cover 
1.1%  of  the  area  in  the  western  United  States  and  are  third  most  extensive  anthropogenic  feature 
in  that  region  (following  agriculture  [9.8%]  and  populated  areas  [1.9%];  Leu  et  al.  2008)  and 
demonstrably  affect  invasion  by  non-native  plants  and  the  number  of  point  sources  of  fires 
(Gelbard  and  Belnap  2003). 


Hierarchical  estimates  of  species  richness 

We  originally  aimed  to  assess  two  methods  for  correcting  bias  in  estimates  of  species  richness  at 
small  spatial  extents:  standardizing  sampling  effort  by  area  sampled  and  by  number  of 
individuals  sampled.  The  relation  between  species  richness  and  area  often  has  one  of  two 
functional  forms:  a  power  function  (linear  on  a  plot  of  log[species  richness]  against  log[area])  or 
an  exponential  function  (linear  on  a  plot  of  species  richness  against  log[area]).  At  least  five 
biological  mechanisms  have  been  hypothesized  to  explain  relations  between  species  richness  and 
area.  The  habitat  diversity  hypothesis  suggests  that  as  area  increases,  the  number  of  niches 
increases.  The  area  hypothesis  is  based  on  the  fact  that  probability  of  extinction  generally 


6 


decreases  as  abundance  increases,  and  abundance  generally  increases  as  area  increases.  The 
passive  sampling  hypothesis  assumes  that  immigration  of  a  given  species  increases  as  area 
increases.  The  resource  concentration  hypothesis  suggests  that  resource  density  increases  as  area 
increases.  The  edge  effects  hypothesis  notes  that  the  perimeter-to-area  ratio  typically  decreases 
as  area  increases. 

Most  published  species  richness-area  curves  were  based  on  one  of  four  sampling  designs  (Figure 
1).  Nested  quadrats  accumulate  area  by  embedding  each  quadrat  in  a  quadrat  of  equal 
proportions  but  larger  size.  Grids  are  contiguous  sample  units  of  the  same  size  and  shape;  area  is 
accumulated  by  increasing  the  number  of  cells  sampled.  Non-contiguous  quadrats  are  a  regular 
array  of  sample  units  of  the  same  size  and  shape,  typically  homogenous  and  embedded  in  a 
homogenous  matrix.  Island-like  patches  are  of  different  size  and  shape,  typically  embedded  in  a 
homogenous  matrix. 


(a)  (b) 
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Figure  1.  Sampling  designs  for  species  richness-area  curves:  nested  quadrats  (a),  grid  (b),  non¬ 
contiguous  quadrats  (c),  island- like  habitat  patches  (d). 
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Species  richness-area  plots  based  on  quadrat  sampling  designs  may  reach  the  same  asymptote, 
but  the  number  of  quadrats  necessary  to  reach  the  asymptote  decreases  as  quadrat  size  increases. 
Thus,  species-area  relations  are  affected  by  both  the  spatial  extent  and  resolution  of  sampling. 

As  we  pursued  this  objective,  we  realized  that  some  of  our  sampling  designs  differed  from  those 
typically  represented  in  the  literature,  which  made  it  difficult  to  apply  existing  species  richness- 
area  models.  We  sampled  anurans  in  island- like  patches  in  a  heterogeneous  matrix.  However,  the 
regional  species  pool  of  anurans  is  relatively  well-known,  and  rarefaction  resulted  in  estimates  of 
species  richness  that  considerably  exceeded  that  of  the  species  pool.  Therefore,  it  did  not  make 
sense  to  apply  species  richness-area  models  to  anurans.  We  sampled  birds  in  all  ecosystems  in 
heterogeneous  patches  of  equal  size  and  shape  that  were  embedded  in  a  heterogeneous  matrix. 
We  sampled  birds  in  this  manner  because  point-count  surveys  are  the  most  common  method  of 
sampling  birds  (see  next  section);  conducting  point-count  surveys  will  allow  our  data  to  be 
incorporated  into  numerous  other  analyses  that  also  are  based  on  point  counts.  Additionally, 
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because  land  cover  and  topography  in  the  ecosystems  in  which  we  worked  is  highly 
heterogeneous,  it  was  not  possible  to  establish  a  sampling  array  of  homogenous  patches  in  a 
homogenous  matrix.  We  also  sampled  butterflies  in  the  Chesapeake  Bay  Lowlands  in 
heterogeneous  patches  of  equal  size  and  shape  that  were  embedded  in  a  heterogeneous  matrix. 
We  sampled  butterflies  in  the  central  and  western  Great  Basin  in  heterogeneous  patches  of 
unequal  size  embedded  in  a  heterogeneous  matrix.  Our  sampling  design  for  butterflies  in  the 
central  and  western  Great  Basin  allowed  us  to  pool  data  from  the  current  project  with  data 
collected  in  previous  years  (this  was  not  a  consideration  for  collection  of  butterfly  data  in  the 
Chesapeake  Bay  Lowlands).  The  latter  sampling  design  also  was  developed  to  accommodate 
extensive  elevational  gradients  in  the  Great  Basin,  which  our  previous  work  suggested  were 
associated  with  species  richness  and  species  occurrence. 

Therefore,  we  developed  both  novel  data-aggregation  methods  and  novel  statistical  methods  to 
address  species  richness-area  relations  in  our  data.  These  methods  also  allowed  us  to  explore 
relations  between  species  richness  and  environmental  gradients  other  than  area,  such  as  elevation 
or  topographic  heterogeneity,  that  may  explain  considerable  variation  in  species  richness.  The 
ability  to  address  other  environmental  gradients  is  highly  relevant  given  landscape  heterogeneity; 
that  is,  other  environmental  attributes  cannot  be  held  constant  while  area  varies. 


Periodicity  of  sampling:  duration  of  point  counts  of  birds 

Our  second  objective  was  to  develop  guidelines  for  periodicity  of  sampling.  Point-count  surveys 
are  the  most  common  method  of  sampling  birds  (Rosenstock  et  al.  2002,  Bart  2005,  Buckland 
2006),  and  many  investigators  have  examined  whether  species-  or  community-level  inferences 
derived  from  point  counts  are  affected  by  sampling  design  (e.g.,  Ralph  et  al.  1993,  1995,  Petit  et 
al.  1995,  Smith  et  al.  1995,  Thompson  and  Schwalback  1995,  Bibby  et  al.  2000,  Shiu  and  Lee 
2003,  Esquivel  and  Peris  2008,  Cimprich  2009,  Reidy  et  al.  2011).  The  temporal  window  for 
effectively  sampling  birds  on  a  given  day  can  be  relatively  short,  sometimes  less  than  four  hours, 
and  the  breeding  season  in  a  given  region  is  often  limited  to  five  to  eight  weeks.  Therefore,  there 
are  trade-offs  among  the  number  of  point-count  locations  (hereafter  points;  Buskirk  and 
McDonald  1995)  sampled,  the  geographic  extent  of  sampling,  and  the  duration  of  point-count 
surveys  (hereafter  counts;  Buskirk  and  McDonald  1995,  Ralph  et  al.  1995,  Vergara  et  al.  2010). 

Estimates  of  detection  probability  and  abundance  generally  increase  as  count  duration  increases, 
but  there  is  considerable  variation  among  species.  For  example,  detection  probabilities  of  14 
songbird  species  that  breed  in  deciduous  forests  in  the  eastern  United  States  increased  from  >  0.4 
to  >  0.9  when  count  duration  increased  from  5  min  to  20  min.  However,  detection  probabilities 
of  57%  of  these  species  varied  among  years  (Dawson  et  al.  1995).  Similarly,  detection 
probability  increased  as  count  duration  increased  from  3  to  6  to  10  min  for  six  songbird  species 
that  breed  in  deciduous  forests  in  the  eastern  United  States  (Buskirk  and  McDonald  1995).  By 
contrast,  mean  abundance  derived  from  10-min  counts  was  higher  than  that  derived  from  6-min 
counts  for  15%  of  13  species  (Thompson  III  et  al.  2002).  Other  studies  found  that  the  number  of 
individuals  increased  as  count  duration  increased  regardless  of  the  period  during  the  morning  in 
which  3-,  6-,  or  10-min  counts  were  conducted  (Buskirk  and  McDonald  1995)  or  the  time  of  year 
in  which  5-,  10-,  15-,  or  20-min  counts  were  conducted  (Smith  et  al.  1998).  In  southwestern 
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France,  the  abundance  of  90%  of  21  breeding  species  was  greater  when  based  on  10-min  than  5- 
min  counts;  the  abundance  of  all  21  species  was  greater  when  based  on  20-min  than  15 -min 
counts  (Bonthoux  and  Balent  2012).  Estimates  of  the  density  of  tropical  species  that  were 
adjusted  for  imperfect  detection  were  13%  greater  when  based  on  10-min  than  2-min  counts  (Lee 
and  Marsden  2008). 

Count  models  suggest  that  the  effect  of  count  duration  on  precision  is  equivocal.  Precision  on  the 
basis  of  the  coefficient  of  variation  (CV)  was  homogenous  among  6-,  8-,  and  10-min  counts 
(Thompson  et  al.  2002),  and  precision  on  the  basis  of  the  standard  error  (SE)  was  homogenous 
among  3-,  6-,  and  10-min  counts  (Buskirk  and  McDonald  1995).  By  contrast,  Smith  et  al.  (1998) 
found  that  SE-based  precision  decreased  as  count  duration  increased  in  5-min  increments  from  5 
to  20  min.  Although  many  studies  investigated  the  effects  of  increasing  count  duration  on 
inferences  from  models  based  on  counts,  it  is  unclear  whether  count  duration  affects  inferences 
from  occupancy  models  (MacKenzie  et  al.  2002). 

Occupancy  models  (MacKenzie  et  al.  2002)  frequently  are  implemented  to  analyze  point-count 
data  (e.g.,  Betts  et  al.  2008,  Saracco  et  al.  2011,  Frey  et  al.  2012).  Typically,  occupancy  models 
use  data  from  repeated  surveys  at  multiple  locations  to  infer  the  proportion  of  locations  or  area 
that  is  occupied  by  a  given  species  and,  optionally,  the  environmental  attributes  associated  with 
occupancy.  Assessments  of  trade-offs  between  the  number  of  sampling  locations  and  the  number 
of  surveys  at  each  location  have  resulted  in  guidelines  for  sampling  design  (MacKenzie  et  al. 
2002,  MacKenzie  and  Royle  2005,  MacKenzie  et  al.  2006,  Bailey  et  al.  2007).  The  effect  of 
count  duration  on  inferences  from  occupancy  models,  however,  has  received  relatively  little 
attention. 

We  explored  the  effects  of  5-min  and  8-min  count  durations  on  inferences  about  detection 
probability  and  occupancy,  and  precision  of  occupancy  estimates,  during  the  breeding  seasons  of 
2012  and  2013  in  the  Chesapeake  Bay  Lowlands,  central  Great  Basin,  and  western  Great  Basin. 
To  extend  the  comparison  beyond  two  years,  we  also  examined  whether  estimates  of  occupancy 
differed  from  2009  through  2013  for  three  species  in  the  central  Great  Basin  (data  collection 
from  2009  through  2011  was  supported  by  other  sources).  Furthermore,  we  determined  the 
percentage  of  species  for  which  occupancy  models  converged,  i.e.,  diagnostics  indicated  no 
problems  with  parameter  estimation,  and  for  which  detection  probabilities  derived  from  both  5- 
min  and  8-min  counts  were  >  0.3.  The  0.3  threshold  was  suggested  as  a  reasonable  means  to 
minimize  bias  (i.e.,  deviation  of  estimated  occupancy  from  the  true  occupancy  value)  in 
occupancy  estimates  (MacKenzie  et  al.  2002).  We  chose  the  5-min  count  duration  to  be 
consistent  with  standard  point-count  protocols  (Ralph  et  al.  1993,  1995;  Matsuko  et  al.  2014). 

We  chose  the  8-min  duration  because  a  pilot  survey  in  the  Chesapeake  Bay  Lowlands  revealed 
that  about  94%  of  species  were  detected  within  8  min  and  because  8  min  was  the  maximum 
count  duration  at  which  no  reduction  in  the  number  of  points  was  necessary  given  logistic 
constraints  and  associated  travel  time  between  points  in  the  central  and  western  Great  Basin.  The 
geographic  and  temporal  breadth  of  our  analysis  is  novel,  and  our  inferences  are  relevant  to  the 
design  of  surveys  based  on  counts  and  estimates  of  species  richness  based  on  occupancy  models 
(e.g.,  Iknayan  et  al.  2014). 
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Periodicity  of  sampling:  trade-offs  between  single-day  and  multiple-days  sampling 

Many  sampling  and  analytical  methods  have  been  developed  to  address  imperfect  detection  (e.g., 
Buckland  et  al.  2001,  Williams  et  al.  2002).  The  lag  time  between  surveys  is  a  key  consideration 
in  the  design  of  occupancy  studies.  Generally,  one  either  visits  each  site  on  multiple  days  and 
conducts  a  single  survey  on  each  visit  (hereafter,  the  multiple-days  design),  or  visits  each  site  on 
one  day  and  conducts  multiple  surveys  on  that  day  (hereafter,  the  single-day  design;  MacKenzie 
et  al.  2006).  Both  the  multiple-days  and  single-day  designs  have  been  used  to  sample  animals  in 
diverse  ecosystems,  although  the  multiple-days  design  is  more  common.  The  single-day  design, 
for  example,  was  used  to  develop  species  distribution  models  for  birds  that  breed  in  forests  of  the 
eastern  United  States  (Schwenk  and  Donovan  2011).  The  multiple-days  design  was  used  to 
examine  mammalian  responses  to  roads  (Nicholson  and  Van  Manen  2009),  changes  in 
occupancy  of  12  anuran  species  from  2002-2006  (Walls  et  al.  2011),  population  declines  in 
amphibians  from  2002-2011  (Adams  et  al.  2013),  avian  responses  to  vehicular  sound  (Goodwin 
and  Shriver  201 1),  and  detection  probabilities  for  butterflies  (Pellet  2008).  There  are  many 
variations  on  the  single-day  and  multiple-days  designs. 

The  degree  to  which  species  are  falsely  detected  should  not  differ  between  the  single-day  and 
multiple  visits  designs  because  false  positives  reflect  observer  error  (McClintock  et  al.  2010a  and 
2010b).  The  degree  to  which  the  assumptions  of  closure,  heterogeneity  in  detection  probability 
and  occupancy,  and  independence  among  detections  of  a  given  species  are  met,  however,  likely 
differs  between  the  single-day  and  the  multiple-days  designs,  and  these  differences  may  affect 
estimates  of  detection  and  occupancy  probabilities  and  their  precision. 

A  single-day  design  is  more  likely  to  meet  the  closure  assumption  than  the  multiple-days  design 
because  all  surveys  are  completed  over  a  relatively  short  period  of  time  (MacKenzie  et  al.  2006, 
Rota  et  al.  2009).  For  example,  males  of  some  bird  species  vacate  territories  early  in  the  breeding 
season  and  establish  new  territories  for  the  remainder  of  the  breeding  season  (Newton  2000, 

Betts  et  al.  2008).  If  a  species  leaves  a  sample  unit  during  the  period  of  sampling,  then  the 
closure  assumption  is  violated.  Violation  of  the  closure  assumption  will  not  necessarily  cause 
bias  in  estimates  of  occupancy  probability.  If  movements  on  and  off  sample  units  by  the  focal 
species  occur  at  random,  then  bias  in  the  estimator  is  expected  to  be  minor  (MacKenzie  et  al. 
2006,  Kendall  and  White  2009).  However,  if  movements  are  Markovian  (i.e.,  the  presence  of  the 
focal  species  on  the  sample  unit  during  one  survey  is  conditional  on  its  presence  during  previous 
surveys),  then  bias  may  be  greater  (MacKenzie  et  al.  2006,  Kendall  and  White  2009). 

Variation  in  environmental  conditions  (e.g.,  weather)  may  cause  detection  probability  to  differ 
among  surveys  at  a  given  site  or  among  sites  (MacKenzie  et  al.  2006,  Simons  et  al.  2007). 
Although  survey- specific  covariates  can  account  for  some  of  the  heterogeneity,  the  use  of  a 
single-day  design  may  cause  more  heterogeneity  in  detection  probability  among  sites  than  the 
multiple- visit  design.  Consider  a  scenario  in  which  10  sites  are  surveyed  on  the  first  day  of 
sampling  under  a  single-day  design  and  a  set  of  10  different  sites  are  surveyed  on  the  second  day. 
If  weather  conditions  differ  between  the  two  days  in  ways  that  affect  detectability  of  the  species, 
then  detection  probability  at  sites  surveyed  on  the  first  day  will  differ  from  detection  probability 
at  sites  surveyed  on  the  second  day.  The  occupancy  model  is  analogous  to  capture -recapture 
models  for  closed  populations  (MacKenzie  et  al.  2002),  and  heterogeneity  in  capture  probability 
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negatively  biases  estimators  of  abundance  from  capture-recapture  models  (Kendall  1999). 
Similarly,  we  expect  that  if  the  single-day  design  induces  heterogeneity  in  detection  probability 
among  sites,  then  estimators  of  occupancy  probability  will  be  biased  negatively  (MacKenzie  et 
al.  2006).  When  the  multiple-days  design  is  implemented,  weather  conditions  and  any  associated 
detection  biases  are  likely  to  differ  both  among  sites  and  among  surveys  at  a  given  site.  Because 
all  sites  in  the  sample  rarely  can  be  surveyed  on  the  same  day,  weather  conditions  are  not 
replicated  among  sites.  Nevertheless,  conventional  wisdom  suggests  that  the  multiple-days 
design  is  more  likely  than  the  single-day  design  to  capture  the  range  of  potential  weather 
conditions  within  the  bounds  of  sampling  constraints  to  which  the  population  of  sites  is  exposed 
during  the  season. 

The  multiple-days  design  is  more  likely  to  meet  the  assumption  of  independent  observations  than 
the  single-day  design  because  surveys  typically  are  separated  by  a  longer  period  of  time.  An 
observer  is  more  likely  to  remember  species  detected  during  previous  surveys  at  a  site  if 
relatively  little  time  elapses  between  surveys,  particularly  if  the  same  observer  conducts  all 
surveys.  Whether  surveys  are  independent  can  be  tested  by  including  information  on  previous 
detections  (i.e.,  capture  history  or  detection  history)  as  a  covariate  in  a  model  of  detection 
probability,  and  examining  whether  detection  probability  was  higher  at  sites  where  the  species 
previously  was  detected  (MacKenzie  et  al.  2006). 

We  used  occupancy  data  on  birds  and  anurans  in  the  Chesapeake  Bay  Lowlands  and  butterflies 
in  the  Great  Basin  to  compare  estimates  of  occupancy  probabilities  that  were  based  on  data  from 
the  two  designs.  To  evaluate  potential  mechanisms  for  differences  in  the  estimates,  we  also 
examined  estimates  of  detection  probability,  detection  histories,  and  support  for  the  detection- 
history  model  from  either  design. 


Estimation  of  species  richness  of  multiple  taxonomic  groups  on  the  basis  of  a  single  group 

As  reflected  in  the  statement  of  need  to  which  our  work  responded,  interest  in  identifying 
standard  and  affordable  ways  to  estimate  species  richness — the  number  of  species  in  a  given 
location  and  time  period — has  been  maintained  for  decades.  Methods  for  estimating  species 
richness  on  the  basis  of  incidence  or  abundance  data  include  rarefaction  (Sanders  1968)  and  the 
Chao  1  (Chao  1984),  Chao  2  (Chao  1987),  and  jackknife  (Burham  and  Overton  1978,  1979) 
estimators.  Several  more-recent  methods  of  estimating  species  richness  capitalize  on  increases  in 
computing  power  or  the  popularity  of  occupancy  models.  For  example,  stacked  and  joint  species 
distribution  models  (SDMs)  (Dubuis  et  al.  2011,  Pollock  et  al.  2014)  make  use  of  single-species 
SDMs,  which  project  probability  of  occurrence  on  the  basis  of  environmental  variables  that  are 
measured  in  the  field  or  derived  from  remotely  sensed  data  (Elith  and  Leathwick  2009). 
Multiple-season,  multiple-species  hierarchical  Bayesian  models  (Gelman  and  Hill  2007,  Dorazio 
et  al.  2010)  are  another  example  of  the  many  relatively  new,  computationally  intensive  methods 
that  can  be  used  to  estimate  species  richness  and  to  explain  variation  in  occupancy  as  functions 
of  environmental  co variates,  ideally  predicting  occupancy  in  space  or  time. 

In  contrast  to  these  newer  methods,  over  time,  many  researchers  have  attempted  to  identify  a 
small  set  of  species  with  occurrence  patterns  that  are  related  to  species  richness  of  a  larger  set  of 
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organisms — i.e.,  indicator  species  (Pearson  1994,  Morrison  et  al.  2012).  Others  have  defined 
indicator  species  as  species  that  are  characteristic  of  land-cover  types  or  environmental 
conditions  (e.g.,  Niemi  and  McDonald  2004,  De  Caceres  et  al.  2010).  The  latter  definition,  and 
associated  methods  (e.g.,  Dufrene  and  Legendre  1997),  differ  from  those  we  use  here. 

If  external  evaluations  suggest  that  the  indicator  species-based  models  accurately  predict  species 
richness,  then  it  may  be  possible  to  monitor  the  occurrence  of  a  small  number  of  indicator 
species  rather  than  conducting  comprehensive  species  inventories.  In  many  temperate 
ecosystems,  the  cost  of  sampling  the  occurrence  of  all  species  of  certain  taxonomic  groups,  such 
as  passerines,  is  not  appreciably  greater  than  that  of  sampling  a  small  subset  of  the  group.  For 
taxonomic  groups  that  are  more  cryptic,  monitoring  a  subset  of  an  assemblage  may  be  much 
more  feasible  than  monitoring  the  full  assemblage.  Identification  of  indicator  species  can  be 
simpler  computationally  than  innovative  but  largely  unproven  methods  such  as  the 
implementation  of  SDMs  or  hierarchical  Bayesian  models,  and  the  relations  between  indicator 
species  and  species  richness  may  be  easier  to  interpret  and  to  communicate  to  end-users. 

We  conducted  work  on  indicator  species  in  the  central  Great  Basin  before  the  current  project  was 
initiated.  In  that  earlier  work,  we  used  objective  statistical  methods  to  identify  butterflies  and 
birds  that  could  indicate  species  richness  of  the  same  taxonomic  group  or  the  other  taxonomic 
group  (Mac  Nally  and  Fleishman  2002,  2004,  Thomson  et  al.  2005,  2007).  For  example,  a  model 
based  on  the  occurrence  patterns  of  five  species  of  butterflies  explained  88%  of  the  deviance  of 
species  richness  of  56  butterfly  taxa;  when  predictions  of  this  model  were  confronted  with  new 
data  from  a  nearby  set  of  locations,  more  than  90%  of  the  observed  values  fell  within  the  95% 
credible  intervals  of  the  predictions  (Mac  Nally  and  Fleishman  2004).  In  our  previous  work,  we 
also  built  models  with  data  on  birds  and  butterflies  that  were  collected  from  1996-2003  in  three 
mountain  ranges  in  the  central  Great  Basin,  and  used  bootstrapping  to  conduct  internal 
evaluations  of  the  models.  Furthermore,  in  our  previous  work,  we  used  new  data  on  birds, 
collected  in  2004  in  25  previously  unsampled  locations  in  one  of  the  central  Great  Basin 
mountain  ranges,  to  conduct  a  preliminary  external  evaluation  of  the  models  (Thomson  et  al. 
2007),  but  we  did  not  have  sufficient  data  to  explicitly  test  the  transferability  of  indicator  species. 

The  extent  to  which  particular  indicator  species  are  transferable  is  relevant  to  the  management  of 
extensive  areas,  including  those  in  which  there  may  be  geographic  differentiation  in  the 
responses  of  widely  distributed  organisms  to  their  environment.  Therefore,  in  the  current  project, 
we  sought  to  develop  a  transferable  model  for  identifying  sets  of  indicator  species  that  explained 
considerable  variation  in  species  richness.  In  this  project,  we  also  compared  predictions  of 
species  richness  based  on  these  indicator  species  to  predictions  based  on  a  set  of  environmental 
variables.  Our  current  project  allowed  us  to  conduct  external  evaluations  of  the  spatial  and 
temporal  transferability  of  indicator  species  and  environmental  variables  that  were  considerably 
more  rigorous  than  our  previous  (Thomson  et  al.  2007)  evaluations.  In  the  current  work,  for 
example,  we  tested  the  extent  to  which  indicator  species  identified  from  models  built  with  data 
from  the  central  Great  Basin  predicted  variation  in  species  richness  of  birds  and  butterflies  in  the 
zoogeographically  distinct  western  Great  Basin  (Behle,  1963,  1978;  Austin  and  Murphy  1987), 
and  vice  versa.  In  our  current  work,  we  also  tested  whether  models  built  from  data  collected  over 
one  to  three  years  in  a  given  subregion  explained  variation  in  species  richness  in  that  same 
subregion  during  one  to  three  years  of  later  sampling. 
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Theory  and  practical  application  of  occupancy  models:  hierarchical  models  of  occupancy 

Species  richness  traditionally  has  been  measured  as  the  number  of  species  observed.  However, 
species  richness  also  can  be  represented  as  the  aggregate  probability  of  occurrence  of  many 
individual  species.  Relatively  few  studies  have  combined  single-species  occupancy  models 
to  estimate  patterns  of  occupancy  of  multiple  species  and  to  relate  these  community-  or 
assemblage-level  patterns  to  environmental  variables  (but  see  Carillo-Rubio  et  al.  2014).  In 
theory,  multiple-species  occupancy  models  (Royle  and  Dorazio  2008)  eliminate  the  need  to 
run  many,  potentially  dozens,  of  single-species  models.  Hierarchical  occupancy  models  are 
applicable  to  analysis  of  nested  data  (e.g.,  species  occurrences  nested  within  species  richness 
of  a  community)  (Gelman  and  Hill  2007)  because  they  allow  for  modeling  relations  among 
levels  in  the  hierarchy  and  variability  in  species  richness  at  all  levels.  For  example,  a 
hierarchical  model  allows  one  to  account  for  variation  in  detection  probabilities  among 
species  when  estimating  credible  intervals  around  species-richness  metrics  at  the  community 
level.  Like  all  occupancy  models,  hierarchical  models  also  allow  one  to  differentiate  the 
sampling  process  (detection)  and  the  state  process  (occupancy  of  the  community). 

Hierarchical  models  assume  that  all  species  within  a  community  are  members  of  one 
multiple- species  population.  Therefore,  these  models  allow  one  can  obtain  occupancy 
estimates  for  species  with  detection  probabilities  that  are  too  low  to  enable  single-species 
occupancy  estimation.  In  these  cases,  species  with  relatively  high  detection  probabilities 
inform  occupancy  estimates  for  those  with  relatively  low  detection  probabilities  (Dorazio  et 
al.  2010).  However,  a  relatively  small  number  of  species  can  disproportionately  affect 
inferences  about  covariates’  associations  with  species  richness  (Zipkin  et  al.  2009). 

We  used  hierarchical  models  of  occupancy  to  explore  species  richness  of  the  community  of 
breeding  birds  in  the  Chesapeake  Bay  Lowlands  and  the  central  and  western  Great  Basin  and 
abiotic  and  biotic  environmental  attributes  associated  with  species  richness.  Because  the 
biology  and  environmental  associations  of  the  species  within  each  community  are  highly 
diverse,  we  grouped  species  into  three  guilds  on  the  basis  of  nesting  stratum:  on  the  ground  or 
in  low  shrubs,  in  tall  shrubs  or  trees,  or  in  tree  cavities. 


Theory  and  practical  application  of  occupancy  models:  estimation  of  the  occupancy  of 
butterflies  in  diverse  biogeographic  regions 

The  population  dynamics  of  individual  species  of  butterflies  have  been  examined  via  both  mark- 
recapture  analyses  (e.g..  Brown  and  Ehrlich  1980,  Fleishman  et  al.  2002,  Leidner  and  Haddad 
2011)  and  occupancy  models  (e.g.,  Pellet  2008,  van  Strien  et  al.  2011,  Bried  et  al.  2012,  Roth  et 
al.  2014).  Collection  of  data  on  occupancy  of  many  species  within  a  butterfly  assemblage,  and 
therefore  application  of  occupancy  models  to  many  species  rather  than  to  one  or  a  small  number 
of  species,  is  complicated  by  many  sources  of  variation  in  phenology  (e.g.,  Baughman  et  al. 

1988,  Weiss  et  al.  1988).  Moreover,  the  number  of  generations  per  year  varies  among  and  within 
species,  and  can  be  plastic.  Accordingly,  it  is  quite  difficult  to  gauge,  a  priori,  the  period  in 
which  a  given  species  is  available  for  sampling. 
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Assemblage-level  surveys  of  butterflies  traditionally  addressed  variation  in  phenology  by 
conducting  surveys  every  one  to  three  weeks  across  the  assemblage’s  flight  season.  This  method 
maximizes  the  likelihood  that  at  least  one  or  two  surveys  will  coincide  with  each  species’  flight. 
However,  one  survey  is  insufficient  to  estimate  detection  probability,  and  over  several  weeks, 
butterfly  assemblages  are  not  closed.  If  the  flight  season  of  a  species  can  be  estimated,  then 
multiple  surveys  potentially  can  be  used  to  develop  a  detection  history.  But  the  occupancy  status 
of  sample  units  is  not  constant  across  the  season,  and  changes  in  occupancy  are  not  random. 
Therefore,  detection  estimates  that  treat  multiple  surveys  as  replicate  samples  of  the  same  species 
may  be  negatively  biased,  and  the  resulting  estimates  of  occupancy  positively  biased.  Another 
possibility,  not  mutually  exclusive,  is  to  sample  each  site  repeatedly  on  each  sampling  date  (the 
robust  design;  Kendall  et  al.,  1997).  A  third  option  is  to  relax  the  closure  assumption  by 
assuming  that  species  are  available  for  sampling  at  different  times  (Kendall  et  al.  2013). 

We  explored  the  extent  to  which  occupancy  of  butterflies  in  the  Chesapeake  Bay  Lowlands, 
central  Great  Basin,  and  western  Great  Basin  could  be  explained  on  the  basis  of  vegetation, 
topography,  and  other  environmental  attributes.  We  included  survey-specific  covariates  of 
detection  probability,  and  we  modeled  occupancy  and  detection  in  multiple  years  as  a  function  of 
the  same  covariates  to  examine  the  temporal  transferability  of  results.  We  also  examined  whether 
results  for  species  that  occur  in  two  of  the  assemblages  were  geographically  consistent. 
Furthermore,  we  examined  the  degree  to  which  assumptions  of  closure  at  the  single- species  level 
were  met  with  an  assemblage-level  sampling  design. 


Theory  and  practical  application  of  occupancy  models:  environmental  associations  with 
multiple  states  of  anuran  occupancy 

The  occupancy  model  of  MacKenzie  et  al.  (2002,  2006)  does  not  address  abundance  at  each 
sample  unit.  This  limitation  to  understanding  of  the  causes  of  spatial  variation  in  the  abundances 
of  populations  led  to  extensions  of  the  model  that  accommodate  data  on  the  relative  abundance 
and  other  attributes  of  species  at  sample  units  (Royle  and  Nichols  2003,  Nichols  et  al.  2007, 
MacKenzie  et  al.  2009).  One  set  of  extensions  was  developed  for  data  from  surveys  that  record 
multiple  states  of  occupancy  (Nichols  et  al.  2007,  MacKenzie  et  al.  2009),  where  states  include 
abundance  or  density  classes  (MacKenzie  et  al.  2009,  Veran  et  al.  2016).  These  models  have 
been  used  to  evaluate  hypothesized  associations  between  attributes  of  sample  units  and  the 
relative  abundances  of  populations  that  occupy  the  sample  units.  For  example,  Goswami  et  al 
(2014)  recorded  the  density  of  elephant  dung  piles  and  tracks  on  sample  units  and  used  multiple- 
state  occupancy  models  to  evaluate  the  roles  of  protected  areas  and  human  presence  on  the 
intensity  of  land  use  by  Asian  elephants  ( Elephas  maximus). 

Concerns  about  declines  in  populations  of  amphibians  in  the  1980s  and  1990s  (Wake  1991)  led 
to  the  development  of  several  state  and  federal  amphibian  monitoring  programs  (e.g.,  the  North 
American  Amphibian  Monitoring  Program  (Weir  and  Mossman  2005).  In  many  of  these 
programs,  the  calling  intensity  of  breeding  choruses  at  wetlands  is  recorded  as  an  ordered 
categorical  variable  (Weir  and  Mossman  2005).  The  data  from  these  surveys  provide  information 
not  only  on  occupancy,  but  on  relative  abundance  across  sample  units.  We  used  similar  field 
methods  and  data  to  identify  attributes  of  wetlands  and  the  surrounding  terrestrial  environment 
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that  were  associated  with  populations  of  high  relative  abundance.  We  compared  those  attributes 
to  the  attributes  of  wetlands  occupied  by  amphibian  populations  of  any  abundance  class. 


Responses  of  songbird  occupancy  to  environmental  change 

Increases  in  human  population  size  and  associated  urbanization  are  a  major  source  of  land-cover 
fragmentation  and  contemporary  environmental  change.  These  changes  can  increase  the  amount 
and  quality  of  habitat  for  some  native  species  of  animals,  including  white-tailed  deer.  The  deer, 
in  turn,  have  the  potential  to  further  change  levels  of  land-cover  fragmentation  and  vegetation 
structure  by  consuming  herbaceous  and  woody  vegetation.  Populations  of  both  native  and  non¬ 
native  deer  have  increased  in  many  locations  worldwide,  prompting  studies  of  ecosystem 
responses  to  abundant  deer  in  Europe  and  North  America  (Knox  1997,  Cote  et  al.  2004). 

Intensive  browsing  by  deer  can  lead  to  simplification  and  reduction  of  biomass  in  understory 
vegetation  (Fuller  2001,  Stockton  et  al.  2005,  Martin  et  al.  2010),  especially  when  deer  become 
less  selective  as  their  abundance  increases  (Augustine  and  McNaughton  1998).  Such  vegetation 
changes  are  believed  to  have  prompted  declines  in  foliage-dwelling  invertebrates  (Allombert  et 
al.  2005).  Vegetation  volume  affects  the  quality  of  habitat  for  many  breeding  songbird  species 
(Mills  et  al.  1991).  Therefore,  high  abundances  of  deer  can  change  habitat  structure  and  food 
availability  for  birds,  especially  those  that  are  associated  with  lower  forest  strata. 

In  the  eastern  United  States,  the  abundance  of  white-tailed  deer  has  increased  substantially  over 
the  past  century  (Cote  et  al.  2004).  For  example,  between  1931  and  early  1990,  the  deer 
population  in  Virginia  was  estimated  to  have  grown  36-fold  as  predators  were  extirpated,  large 
areas  in  which  deer  formerly  were  hunted  became  urbanized,  and  land  use  fragmented  forest 
cover  (McShea  et  al.  2003,  Fovely  et  al.  2013).  Fragmentation  provides  deer  with  ready  access  to 
abundant  forage  at  sunlit  edges  between  forest  and  roadsides,  grasslands,  agricultural  fields,  or 
residential  neighborhoods,  while  providing  cover  in  nearby  forest  (Augustine  and  de  Calesta 
2003,  Fovely  et  al.  2013). 

Several  studies  have  suggested  a  link  between  deer  browsing  and  songbird  dynamics.  Species 
richness  and  abundance  declined  as  deer  density  increased  in  Pennsylvania  (deCalesta  1994)  and 
Massachusetts  (DeGraaf  et  al.  1991).  Islands  in  southern  British  Columbia  with  high  densities  of 
black-tailed  deer  ( Odocoileus  hemionus)  had  lower  densities  of  Rufous  Hummingbirds 
(Selasphorus  rufus).  Song  Sparrows  ( Melospiza  melodia),  and  Fox  Sparrows  ( Passerella  iliaca) 
than  deer-free  islands  (Martin  et  al.  201 1).  Further  north  in  British  Columbia,  on  the  Haida  Gwaii 
archipelago,  long-term  declines  in  abundance  and  diversity  of  understory  birds  were  connected 
to  the  depletion  of  understory  vegetation  by  introduced  black-tailed  deer  (Chollet  et  al.  2014). 
Following  experimental  deer  exclosure  in  Virginia,  abundance  of  ground-nesting  and 
intermediate  canopy-nesting  birds  increased  as  understory  vegetation  regenerated  (McShea  and 
Rappole  2000).  In  Britain,  Eurasian  Blackcaps  ( Sylvia  atricapilla),  which  are  associated  with 
early  successional  forest,  were  more  abundant,  settled  earlier  each  spring,  and  had  better  body 
condition  in  areas  where  roe  deer  ( Capreolus  capreolus)  and  introduced  Reeves’s  muntjac 
( Muntiacus  reevesi )  were  experimentally  excluded  (Holt  et  al.  2013).  At  the  state  and  provincial 
levels,  the  abundance  of  deer  was  related  to  declines  of  understory-nesting  and  foraging  birds 
(Chollet  and  Martin  2013). 
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These  effects  of  deer  on  birds,  primarily  mediated  by  browsing  on  understory  foliage,  can  be 
magnified  in  forest  fragments  if  deer  activity  becomes  concentrated  in  relatively  small  forest 
patches  (Augustine  and  de  Calesta  2003).  Evidence  of  the  effects  of  deer  on  songbirds  primarily 
is  derived  from  studies  on  isolated  or  small  forested  sites.  Only  two  studies  examined  the  effects 
of  deer  on  avian  communities  at  extents  larger  than  forest  patches,  and  we  are  not  aware  of  any 
research  that  examined  associations  between  deer  and  forest  songbirds  at  a  regional  level  with 
fine-grained  estimates  of  the  density  of  deer  and  birds.  By  quantifying  deer  pellet  and  songbird 
densities  across  a  large  region  in  coastal  Virginia  over  four  years,  we  tested  whether  densities  of 
species  that  nest  and  forage  in  the  lowest  stratum  of  forest  foliage  (i.e.,  vegetation  on  the  ground 
and  in  the  understory)  was  negatively  correlated  with  intensity  of  deer  use.  We  also  explored 
whether  relations  between  deer  and  songbirds  were  consistent  along  a  gradient  of  intensity  of 
deer  use.  Our  work  allowed  us  to  examine  how  songbird  densities  may  change  as  the  level  of 
forest  fragmentation  and  relative  intensity  of  deer  use  changes. 


Partitioning  drivers  of  beta  diversity 

Understanding  the  causes  of  variation  in  species  composition  through  space  and  time  (i.e.,  beta 
diversity),  and  relative  variation  in  species  richness  (alpha  diversity)  and  composition,  is  highly 
relevant  to  management.  However,  variation  in  beta  diversity  is  poorly  understood.  For  example, 
differences  in  temporal  variation  in  species  composition  might  inform  the  duration  of  or  the 
inferences  derived  from  monitoring. 

Beta  diversity  can  be  partitioned  into  two  components,  turnover  and  nestedness,  which  reflect 
different  ecological  processes  (Baselga  2010)  (Figure  2).  Turnover  in  space  or  time  occurs  when 
a  species  that  previously  was  present  becomes  absent,  and  a  different  species  becomes  present. 
Nestedness  occurs  when  the  species  present  in  locations  with  relatively  low  species  richness  are 
statistical  subsets  of  the  species  present  in  locations  with  relatively  high  species  richness.  Beta 
diversity  can  be  partitioned  additively  into  turnover  and  nestedness  components,  with  the  sum  of 
turnover  and  nestedness  equal  to  total  beta  diversity  (Baselga  2010).  Baselga  (2010)  defined  the 
Sprensen  dissimilarity  index  as  total  beta  diversity,  the  Simpson  dissimilarity  index  as  beta 
diversity  due  to  turnover,  and  the  difference  between  these  two  indices  as  beta  diversity  due  to 
nestedness. 

Both  turnover  and  nestedness  can  be  defined  spatially  and  temporally.  Species  composition  can 
differ  among  locations  within  a  region,  resulting  in  high  regional  species  richness  even  if 
individual  locations  are  occupied  by  few  species.  Species  composition  can  change  at  one  location 
through  time,  resulting  in  high  species  richness  over  years  or  decades  even  if  that  location  is 
occupied  by  few  species  over  days  or  months  and  species  richness  at  a  given  point  in  time  is 
fairly  consistent.  Both  spatial  and  temporal  beta  diversity  contribute  to  species  richness,  but  few 
studies  have  considered  the  relative  contributions  of  each  to  local  or  regional  species  richness. 
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Temporal  turnover 


Figure  2.  Examples  of  temporal  and  spatial  turnover  and  nestedness.  Different  letters  denote 
different  species;  small  circles  are  points,  transects,  or  canyons;  and  large  bounding  regions 
(thick  black  lines)  are  canyons  or  mountain  ranges.  Temporal  turnover  occurs  when  species’ 
identities  change  through  time.  Temporal  nestedness  occurs  when  species  are  lost  or  gained 
through  time.  Spatial  turnover  and  nestedness  are  defined  similarly,  but  changes  in  species’ 
identities  or  losses  and  gains  of  species  occur  through  space  rather  than  through  time. 


Comparisons  of  beta  diversity  among  spatial  resolutions  can  identify  whether  there  are  consistent 
patterns  in  species  composition  at  any  of  them,  which  might  provide  insight  into  the  natural  and 
anthropogenic  processes  that  affect  the  composition  of  a  given  assemblage.  Such  understanding 
can  inform  selection  of  spatial  extents  and  resolutions  for  management,  including  monitoring. 

For  example,  high  turnover  of  species  among  locations  in  a  region  might  indicate  that  species 
can  move  easily  among  locations  or  have  relatively  plastic  resource  requirements.  By  contrast, 
low  turnover  of  species  among  locations  in  a  region  might  indicate  an  association  of  these 
species  with  fine-resolution  environmental  characteristics. 

Although  it  is  well  established  that  sampling  is  affected  by  environmental  conditions,  and  that 
local  environmental  heterogeneity  (e.g.,  edges  between  land-cover  types)  may  be  associated  with 
beta  diversity,  relatively  little  work  has  considered  whether  and  how  beta  diversity  changes  along 
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environmental  gradients.  For  example,  it  is  unclear  whether  local  climate,  topography,  or 
vegetation  is  associated  consistently  with  temporal  or  spatial  variation  in  faunal  composition  and, 
if  so,  at  which  spatial  resolutions. 


Understanding  how  species  composition  varies  through  space  and  time  may  be  applicable  to  the 
management  and  monitoring  of  extensive  areas.  We  analyzed  spatial  and  temporal  turnover  and 
nestedness  of  bird  and  butterfly  assemblages  in  the  central  and  western  Great  Basin.  We 
compared  estimates  of  turnover  and  nestedness  at  two  spatial  resolutions  (within  canyons  and 
among  canyons)  to  determine  whether  assemblage  composition  was  more  consistent  at  either 
resolution.  We  are  relating  estimates  of  turnover  and  nestedness  to  environmental  variables  to 
evaluate  whether  these  attribute  are  associated  with  beta  diversity.  We  distinguish  our  analyses 
from  others  that  assessed  correlations  between  beta  diversity  and  environmental  diversity  (e.g., 
Mantel  tests)  because  the  latter  focused  on  variation  in  species  composition  rather  than  variation 
in  beta  diversity. 
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Materials  and  Methods 


Field  sampling 

Study  locations.  Field  sampling  was  integral  to  all  of  our  project’s  objectives  and  tasks.  We 
sampled  anurans,  birds,  and  butterflies  in  the  Chesapeake  Bay  Lowlands.  We  sampled  birds  and 
butterflies  in  the  central  Great  Basin  and  western  Great  Basin.  We  did  not  sample  anurans  in  the 
Great  Basin  because  few  anurans  occur  in  the  region.  More  of  our  analyses  and  results  focus  on 
birds  and  butterflies  than  on  anurans  simply  because  anurans  could  not  be  sampled  in  all  regions. 
However,  all  three  taxonomic  groups  often  are  a  focus  of  assessment,  monitoring,  and 
conservation  planning  and  action.  Therefore,  our  work  not  only  informs  methods  for  evaluting 
the  status  and  potential  effects  of  land  use  on  these  diverse  groups,  but  may  be  directly 
applicable  to  efforts  to  maintain  them. 

In  the  Chesapeake  Bay  Lowlands,  our  study  area  included  the  Virginia  Peninsula  between  Toano 
and  Hampton  (Charles,  City,  Henrico,  James  City,  Newport  News,  Williamsburg,  and  York 
Counties,  Virginia)  and  the  Middle  Peninsula  near  West  Point  (King  and  Queen  and  King 
William  Counties,  Virginia)  (Figure  3).  The  Virginia  Peninsula  includes  of  a  network  of  natural 
wetlands  and  ponds  and  anthropogenic  ponds  that  are  embedded  in  a  matrix  of  remnant  upland 
and  bottomland  hardwood  forest,  agricultural  areas,  and  areas  of  low-  to  high-density  housing. 
Our  sampling  included  Joint  Base  Langley-Eustis.  We  focused  on  Joint  Base  Langley-Eustis 
because  it  was  the  most  accessible  installation  in  the  Chesapeake  Bay  Lowlands.  We  conducted 
field  sampling  with  the  area  recognized  as  Fort  Eustis  prior  to  its  consolidation  with  Langley  Air 
Force  Base  in  2010.  The  primary  mission  of  Fort  Eustis  is  Army  transportation  training,  research 
and  development,  engineering,  and  operations,  including  aviation  and  marine  shipping  activities. 
Access  to  other  installations  in  the  Chesapeake  Bay  Lowlands,  including  Camp  Peary  and  Naval 
Weapons  Station  Yorktown  /  Cheatham  Annex,  is  restricted.  We  requested  but  were  denied 
access  to  Fort  A.P.  Hill. 

Depending  on  hydrology  and  soil  drainage,  canopy  trees  in  deciduous  forest  in  the  Chesapeake 
Bay  Lowlands  may  include  red  maple  (Acer  rubrum),  birch  ( Betula  spp.),  black  walnut  (Juglans 
nigra),  sweetgum  ( Liquidambar  styraciflua),  tulip  popular  ( Liriodendron  tulipifera ),  water  or 
black  tupelo  (Nyssa  aquatica,  N.  sylvatica),  sycamore  (Platanus  occidentalis),  and  oaks 
( Quercus  spp.)  (Weakly  2012).  The  canopy  of  upland  coniferous  forests  is  dominated  by  loblolly 
pine  ( Pinus  taeda).  Sweetgum  also  is  a  dominant  species  in  early  successional  stands  (Monette 
and  Ware  1983,  Weakly  2012).  The  canopy  of  upland  deciduous  forest  is  dominated  by 
American  beech  (Fagus  grandifolia),  oaks,  and  American  holly  ( Ilex  opaca).  Loblolly  pine,  tulip 
popular,  and  sweetgum  also  are  present  (Monette  and  Ware  1983,  Weakly  2012). 
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Figure  3.  Department  of  Defense  Installations  within  the  Chesapeake  Bay  Lowlands  and 
locations  of  other  areas  in  which  we  sampled  anurans,  birds,  and  butterflies. 
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Behle  (1963)  recognized  five  centers  of  differentiation  for  birds  in  the  Great  Basin  (Warner, 
Sierra  Nevada,  western  Great  Basin,  eastern  Great  Basin,  Inyo).  He  defined  a  center  of 
differentiation  as  an  expression  of  the  concept  that  “the  area  so  designated  contains  several 
geographic  races  representing  different  species  and  genera,  that  the  ranges  of  the  several  kinds 
coincide  rather  closely,  and  that  the  characters  of  each  kind  are  best  developed  in  the  particular 
region.  These  features  imply  that  the  different  races  have  originated  in  the  area  either  because  of 
the  modifying  influence  of  some  environmental  factor  or  factors  that  have  had  an  effect  on  the 
several  forms,  or  a  barrier  of  some  sort  has  delimited  the  several  kinds  so  that  they  show  a 
common  distributional  pattern”  (page  1168).  Similarly,  Austin  and  Murphy  (1987)  recognized 
six  centers  of  differentiation  for  butterflies  in  the  Great  Basin  (Warner,  Jarbidge,  Central, 
Toiyabe,  Snake,  and  Inyo).  Our  study  locations  in  the  central  Great  Basin  are  within  the  eastern 
(Behle  1963)  or  Toiyabe  (Austin  and  Murphy  1987)  center  of  differentiation,  and  those  in  the 
western  Great  Basin  fall  within  the  Inyo  center  of  differentiation.  Exploring  similarities  and 
differences  in  ecological  patterns  between  centers  of  differentiation  is  relevant  because  the 
management  plans  of  public  departments  and  agencies  sometimes  assume  that  a  single 
management  prescription  will  be  effective  across  the  entire  Great  Basin.  However,  few  studies 
have  addressed  whether  relations  between  species  richness  and  environmental  variables,  or 
between  individual  species  and  environmental  variables,  are  consistent  among  these  centers. 

Our  central  Great  Basin  study  area  included  much  of  the  adjacent  Shoshone  Mountains  and 
Toiyabe,  Toquima,  and  Monitor  Ranges  (Lander,  Nye,  and  Eureka  Counties,  Nevada).  In  the 
western  Great  Basin,  our  study  area  included  part  of  the  east  slope  of  the  Sierra  Nevada  and  the 
adjacent  Wassuk  Range  and  Sweetwater  Mountains  (Mono  County,  California  and  Mineral, 
Douglas,  and  Lyon  Counties,  Nevada).  Our  study  area  also  encompassed  parts  of  the  Hawthorne 
Army  Depot  and  Marine  Corps  Mountain  Warfare  Training  Center  (Figure  4).  Hawthorne  Army 
Depot  stores  conventional  munitions,  demilitarizes  and  disposes  of  unserviceable,  obsolete,  and 
surplus  munitions;  and  maintains  serviceability  through  inspection  and  renovation  to  ensure 
munitions  readiness.  Although  most  of  the  installation  is  flat,  at  relatively  low  elevation,  and 
sparsely  vegetated,  it  also  includes  Mt.  Grant,  which  is  the  highest  (3440  m)  peak  in  the  Wassuk 
Range  and  Mineral  County,  Nevada.  In  part  because  access  to  Mt.  Grant  long  has  been  restricted, 
and  became  highly  restricted  in  2001,  land  use  is  limited  and  habitat  quality  for  many  taxonomic 
groups,  including  birds  and  butterflies,  is  high.  The  Marine  Corps  Mountain  Warfare  Training 
Center  emphasizes  training  of  combat  forces  to  operate  at  high  elevations,  in  complex  terrain, 
and  in  cold  weather.  Most  training  occurs  on  US  Forest  Service  lands  under  a  special  use  permit. 
Other  DoD  installations  in  the  Great  Basin,  such  as  Naval  Air  Station  Fallon,  Hill  Air  Force 
Range,  Wendover  Range,  Deseret  Test  Center,  and  Dugway  Proving  Ground,  are  at  relatively 
low  elevations,  are  sparsely  vegetated,  and  are  inhabited  by  relatively  few  species  of  songbirds 
and  butterflies.  Nellis  Air  Force  Base  largely  is  located  within  the  Mojave  Desert  rather  than  the 
Great  Basin.  We  requested  but  were  denied  access  to  the  Tonopah  Test  Range. 

Among  the  dominant  land-cover  types  in  the  central  and  western  Great  Basin  are  woodlands 
dominated  by  single-leaf  pinyon  ( Pinus  monophylla)  and  juniper  (, Juniperus  osteosperma,  J. 
occidental is),  shrubsteppe  dominated  by  sagebrush  ( Artemisia  spp.),  and  riparian  woodlands 
dominated  by  deciduous  trees  (e.g.,  aspen  [Populus  tremuloides ],  chokecherry  [Prunus 
virginiana ])  and  shrubs  (e.g.,  willow  [Salix  spp.],  Woods’  rose  [Rosa  woodsii ]).  Jeffrey  pine 
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( Pinus  jeffreyi),  lodgepole  pine  ( Pinus  contorta ),  and  red  fir  ( Abies  magnified )  also  are  dominant 
trees  in  some  parts  of  the  east  slope  of  the  Sierra  Nevada  and  Wassuk  Range. 


NAtonal 


Jarbidge 

Subregion 


Central 

Subregio 


i  mas 


Toiyabe 


Inyb. 

Subregion 


Kilometers 


Mojave 

Subregion 


Twino 

Foils 


I  HIM, tv  .  r 


Warner  ; 
Subregibn 


Foivtt 


Lcracto^L 


Snake 

Subregion 


!StO 


?  Scut  'je>: 


Western  study  sites 
Central  study  sites 
DoD  Installations 


Sources  Em  DeLorme.  NAVTEQ,  TomTom,  Intermap  incre  nent  P  Corp  ,  Ff.] 
GEBCO.  USbs.  FAO^PS/^RS^N,  &&)Base  IGN  Kadas  er  NL.  Ordnance; 
Survey.  Esrl  Jap^a  METl.  Esri  China  (Hong  KorT^lt^isstoRo  ,and  the.GIS  User 
c°mmunity  \  i  V  ' 


Figure  4.  Centers  of  faunal  differentiation  with  the  Great  Basin  (Austin  and  Murphy  1987), 
locations  of  our  study  sites  within  the  Inyo  Subregion  (western  Great  Basin)  and  Toiyabe 
Subregion  (central  Great  Basin),  and  locations  of  most  DoD  installations  in  the  Great  Basin. 
Open  black  circles  indicate  the  locations  of  the  Marine  Corps  Mountain  Warfare  Training  Center 
(left)  and  Mt.  Grant  (right). 
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Anurans.  We  delineated  an  anuran  study  area  of  3768  km2  that  included  the  Virginia  Peninsula 
and  portions  of  the  Richmond  municipal  area.  This  area  was  bounded  by  the  Chesapeake  Bay  to 
the  east,  York  River  to  the  north,  and  James  River  to  the  south.  Approximately  3.4%  of  the  study 
area  was  covered  by  wetlands;  freshwater  emergent  wetlands  covered  1.9%  (US  Fish  and 
Wildlife  Service  2010).  Fifty-four  percent  of  the  study  area  was  classified  as  forest  (SE  GAP 
2010).  Developed  areas  and  agriculture  covered  24.0%  and  12.1%  of  the  study  area,  respectively. 
We  stratified  the  study  area  with  the  smallest  US  Geological  Survey  Hydrological  Unit 
Classification  (HUC  12)  (Seaber  et  al.  1987)  and  2000  census  data  (US  Census  Bureau  2010). 

We  assigned  each  HUC  to  one  of  four  housing-density  strata  (Radeloff  et  al.  2005):  very  low  (> 

0  and  <6.18  houses/km2),  low  (>  6.18  and  <  49.42),  medium  (>  49.42  and  <  741.32),  and  high 
(>  741.32).  We  created  a  sampling  frame  of  wetlands  and  ponds  (hereafter,  wetlands)  by 
combining  data  from  the  National  Wetlands  Inventory  (US  Fish  and  Wildlife  Service  2010)  with 
data  on  a  set  of  wetlands  that  we  had  identified.  Prior  to  the  field  season  in  201 1,  we  used  aerial 
imagery  from  the  2009  Virginia  Base  Mapping  Program  (Virginia  Information  Technologies 
Agency  2009)  to  hand-digitize  wetlands  not  captured  by  the  National  Wetlands  Inventory  in  a 
Geographic  Information  System  (GIS)  (ArcMap  9.0,  ESRI  2009).  From  the  sampling  frame,  we 
randomly  selected  wetlands  from  each  of  the  housing-density  strata;  selected  wetlands  were  >1 
km  apart.  The  accessibility  of  the  wetlands  was  evaluated  prior  to  the  field  season.  If  a  given 
wetland  was  not  accessible  (e.g.,  unsafe  to  access  or  on  private  property  with  restricted  access), 
we  randomly  selected  another  wetland.  The  mean  distance  between  neighboring  wetlands  in  the 
final  sample  was  2.63  km  (standard  deviation  [SD]  1.43). 

Our  sampling  protocol  for  anurans  was  based  on  the  North  American  Amphibian  Monitoring 
Protocol  (Weir  and  Mossman  2005).  We  conducted  surveys  from  February  through  July  to 
capture  the  three  periods  during  which  different  species  breed  (Garrett  2002).  Each  survey 
focused  on  the  peak  month  of  vocalization  for  a  breeding  period  (mid  February  through  mid 
March,  mid  April  through  mid  May,  and  mid  June  through  mid  July).  We  adjusted  the  start  date 
of  each  survey  on  the  basis  of  annual  spring  temperatures  and  precipitation.  We  sampled 
adjacent  to  a  pond  or  wetland.  We  sampled  each  pond  three  times  during  each  breeding  period, 
for  a  total  of  nine  visits  per  pond  per  year.  Three  visits  per  breeding  season  is  standard  for  anuran 
surveys  (Weir  et  al.  2005,  2009).  Visits  began  30  min  after  sunset  and  continued  until  0100 
(Bowers  et  al.  1998,  Bridges  and  Dorcas  2000).  We  noted  the  time  that  each  visit  began.  During 
each  visit,  we  noted  time  of  first  detection  for  all  species  calling  within  5  min  [detection 
probabilities  increase  slightly  as  duration  of  sampling  increases  (Dorcas  et  al.  2009)],  estimated 
abundance  of  males  categorically  (1,  individuals  can  be  counted;  2,  calls  of  individuals  can  be 
distinguished  but  there  is  some  overlap  of  calls;  3,  full  chorus  with  calls  that  are  constant, 
continuous,  and  overlapping),  and  estimated  covariates  known  to  affect  detection  of  anurans 
(e.g.,  ambient  sound,  temperature,  wind  speed). 

At  the  time  of  each  survey,  which  we  recorded,  we  measured  ambient  temperature  (°C)  and  wind 
speed  (m  sec'1)  with  a  hand-held  weather  station  (Kestrel  2000,  Nielsen-Kellerman  Co., 
Boothwyn,  Pennsylvania)  and  wind  speed  in  the  canopy  according  to  the  Beaufort  scale  (see 
www.spc.noaa.gov/faq/tornado/beaufort.html).  We  also  recorded  traffic  volumes  during  the 
count,  which  affect  ambient  sound  levels  and  detection  probabilities  (0,  none;  1,  1  car  passing;  2, 
2-5  cars  passing;  3,  continuous  traffic  passing  [about  6-10  cars];  4,  continuous  traffic  passing, 
with  construction  or  industrial  sounds).  We  did  not  conduct  surveys  if  wind  speed  exceeded  3  on 
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the  Beaufort  scale  (leaves  and  small  twigs  constantly  moving,  light  flags  extended),  during  snow, 
during  rainfall  of  sufficient  intensity  to  affect  hearing  ability,  or  if  temperatures  were  before  the 
threshold  defined  by  the  North  American  Amphibian  Monitoring  Protocol  (5.6°  during  the  first 
breeding  period,  12.8°  during  the  second  breeding  period,  and  18.3°  during  the  third  breeding 
period). 

Anuran  data  were  included  in  analyses  of  periodicity  of  sampling  (trade-offs  between  single-day 
and  multiple-days  sampling)  and  of  environmental  associations  with  multiple  states  of 
occupancy. 

Birds.  To  survey  birds  in  the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin, 
we  conducted  point  counts  ( counts )  during  the  breeding  season  (late  May  through  June)  (Dobkin 
and  Rich  1998).  Despite  differences  in  elevation  and  latitude,  the  timing  of  the  breeding  season 
is  similar  among  locations  and  years.  The  sampling  period  overlaps  among  ecosystems  in  part 
because  of  the  protracted  songbird  migration  in  the  eastern  United  States.  The  peak  of  the 
migration  period  in  the  eastern  United  States  varies  among  years  but  generally  occurs  during  the 
first  two  weeks  of  May.  We  did  not  survey  birds  in  the  Chesapeake  Bay  Lowlands  until  most 
species  that  breed  further  north  had  passed  through  the  area.  Similarly,  birds  that  breed  in  the 
Great  Basin  are  present  in  that  region  earlier  in  May,  but  the  ratio  of  birds  that  breed  in  the  Great 
Basin  to  birds  that  breed  further  north  increases  throughout  the  month. 

Most  point  centers  were  >  350  m  apart,  which  minimized  the  probability  that  observer  activity 
would  cause  individuals  to  move  among  points.  During  each  visit,  we  recorded  by  sound  or  sight 
all  birds  using  terrestrial  habitat  within  the  point.  We  visited  each  point  three  times  per  year  for  8 
min  per  count  (Buckland  et  al.  2001,  Siegel  et  al.  2001,  Dickson  et  al.  2009).  We  recorded 
whether  a  detection  occurred  during  the  first  5  min  or  the  last  3  min.  In  each  of  our  ecosystems, 
observers  used  a  laser  rangefinder  to  estimate  the  distance  to  each  bird  detected.  In  the 
Chesapeake  Bay  Lowlands,  we  estimated  distance  as  a  continuous  variable.  In  the  Great  Basin, 
we  sampled  birds  with  fixed-radius  point  counts,  and  used  distance  bins:  <10  m,  >10-25  m, 
>25-50  m,  >50-75  m,  >75-100  m,  and  >  100  m.  The  distance  bins  and  distance  thresholds  we 
used  followed  standard  protocols  outlined  by  Ralph  et  al.  (1993,  1995)  and  reiterated  by 
Matsuoka  et  al.  (2014).  In  our  analyses,  we  included  detections  within  100  m  of  the  observer,  but 
generally  did  not  differentiate  distances  within  that  radius.  We  maintained  and  archived  records 
of  birds  detected  at  distances  of  >100  m. 

At  the  time  of  each  count  in  the  Chesapeake  Bay  Lowlands  (we  recorded  the  time),  we  recorded 
ambient  temperature  (°C)  and  wind  speed  (m  sec"1)  with  a  hand-held  weather  station  (Kestrel 
2000,  Nielsen-Kellerman  Co.,  Boothwyn,  Pennsylvania)  and  wind  speed  in  the  canopy  according 
to  the  Beaufort  scale.  We  also  recorded  traffic  volumes  during  the  count,  which  affect  ambient 
sound  levels  and  detection  probabilities  (0,  none;  1,  1  car  passing;  2,  2-5  cars  passing;  3, 
continuous  traffic  passing  [about  6-10  cars];  4,  continuous  traffic  passing,  with  construction  or 
industrial  sounds). 

In  the  Chesapeake  Bay  Lowlands,  we  established  points  within  riparian  forests  and  upland 
coniferous  and  deciduous  forests.  In  the  central  and  western  Great  Basin,  points  were  positioned 
to  sample  the  dominant  land-cover  types  throughout  the  canyons. 
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The  common  names  of  birds  that  occur  in  the  United  States  are  standardized  by  the  North 
American  Classification  Committee,  an  official  committee  of  the  American  Ornithological 
Society.  Our  use  of  common  and  scientific  names  of  birds  follows  the  Birds  of  North  and  Middle 
America  Checklist ,  the  official  source  on  the  taxonomy  of  birds  in  North  and  Middle  America, 
including  adjacent  islands. 

Bird  data  from  the  Chesapeake  Bay  Lowlands  were  included  in  analyses  of  periodicity  of 
sampling  (duration  of  point  counts),  hierarchical  occupancy,  and  responses  of  songbird 
occupancy  to  environmental  change. 

Bird  data  from  the  central  and  western  Great  Basin  were  included  in  analyses  of  hierarchical 
species  richness,  periodicity  of  sampling  (duration  of  point  counts  and  trade-offs  between  single¬ 
day  and  multiple-days  sampling),  estimation  of  species  richness  of  multiple  taxonomic  groups  on 
the  basis  of  a  single  group,  hierarchical  occupancy,  and  drivers  of  beta  diversity. 

Butterflies.  We  sampled  butterflies  along  transects,  all  of  which  encompassed  points  at  which 
we  surveyed  birds.  Because  timing  and  duration  of  flight  varies  among  butterfly  species, 
locations,  and  years  (Scott  1986),  we  sampled  each  transect  approximately  every  two  weeks 
throughout  the  majority  of  the  flight  season  (late  May  through  mid  August).  During  every  visit, 
we  walked  the  length  of  each  transect  at  a  near-constant  pace  and  recorded  all  butterfly  species 
detected  (Pollard  and  Yates  1993).  Walking  along  transects  is  a  standard  method  for  surveying 
butterfly  assemblages  (Pollard  and  Yates  1993,  Pullin  1995).  Species  that  cannot  be  identified  on 
the  wing  are  captured  and  either  identified  in  the  hand  or  collected  for  identification  in  the 
laboratory. 

In  the  Chesapeake  Bay  Lowlands,  we  sampled  butterflies  along  a  0.5  km  transect  within  each 
bird  point.  In  the  central  Great  Basin,  we  divided  each  canyon  into  100-m  vertical  elevation 
bands  from  its  mouth  to  its  crest.  In  the  western  Great  Basin,  transects  for  butterflies  covered  the 
length  of  each  canyon  and  were  centered  on  bird  points.  We  estimated  abundance  of  each  species 
along  each  transect.  During  each  visit,  we  recorded  the  relative  abundance  (none,  low,  moderate, 
high)  of  individual  plants  (primarily  forbs)  from  which  one  or  more  species  of  butterflies  in  those 
ecosystems  are  known  to  take  nectar  and,  in  the  Great  Basin,  the  relative  abundance  of  sources  of 
mud,  such  as  stream  crossings  (none,  low,  moderate,  high).  Female  fecundity  in  some  species  is 
related  to  nectar  volume  (Boggs  and  Ross  1993),  and  many  species  feed  on  dissolved  minerals  in 
moist  soil  (Scudder  1889,  Arms  et  al.  1974). 

Nomenclature  for  butterflies  in  the  Chesapeake  Bay  Lowlands  follows  Pelham  2016. 
Nomenclature  for  butterflies  in  the  Great  Basin  follows  Fleishman  et  al.  1997  and  Austin  1998. 

Butterfly  data  from  the  Chesapeake  Bay  Lowlands  were  included  in  analyses  of  the  occupancy  of 
butterflies  in  diverse  biogeographic  regions.  Butterfly  data  from  the  central  and  western  Great 
Basin  were  included  in  analyses  of  hierarchical  species  richness,  periodicity  of  sampling  (trade¬ 
offs  between  single-day  and  multiple-days  sampling),  estimation  of  species  richness  of  multiple 
taxonomic  groups  on  the  basis  of  a  single  group,  occupancy  of  butterflies  in  diverse 
biogeographic  regions,  and  drivers  of  beta  diversity. 
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Vegetation  and  other  environmental  co variates.  To  characterize  vegetation  composition  and 
structure  at  each  pond  where  we  sampled  anurans,  we  established  three  1  x  3  m  sampling  plots 
along  the  edge  of  the  pond  (Figure  5).  One  plot  was  positioned  at  the  location  at  which  anurans 
were  sampled.  The  other  two  plots  were  positioned  at  20  m  on  either  side  of  the  anuran- sampling 
location.  In  each  plot  we  measured  water  pH  and  water  depth  at  1  m  from  the  edge  of  the  pond. 
We  measured  composition  and  structure  of  aquatic  and  terrestrial  vegetation  in  three  adjacent  1- 
m2  quadrats.  One  quadrat  overlapped  the  water-land  edge,  one  extended  into  the  water,  and  one 
extended  onto  the  land.  In  each  plot  we  measured  the  percentage  of  cover  of  trees  >  3  m,  shrubs 
of  different  heights  (>  0.5-3  m,  >  0.3-0. 5  m,  0. 1-0.3  m,  and  <0.1  m),  grasses,  forbs,  moss  and 
lichens,  floating  vegetation,  emergent  vegetation,  and  submerged  vegetation  (Mazerolle  et  al. 
2005).  We  estimated  the  number  of  days  without  precipitation  and  the  number  of  days  since 
above-average  precipitation  (average  within  each  of  the  three  breeding  periods)  on  the  basis  of 
records  from  the  nearest  airport.  We  also  calculated  the  number  of  minutes  between  sunset  and 
initiation  of  the  survey  and  the  Julian  date  (standardized  to  the  first  day  of  the  survey). 


Figure  5.  Design  of  vegetation  sampling  for  anurans.  Figure  not  to  scale. 


Data  on  vegetation  and  other  environmental  co  variates  associated  with  anurans  were 
incorporated  into  analyses  of  environmental  associations  with  multiple  states  of  occupancy. 

Methods  for  sampling  local  vegetation  composition  and  structure  at  bird-survey  points  differed 
between  the  Chesapeake  Bay  Lowlands  and  the  Great  Basin  (we  used  the  same  methods  in  the 
central  and  western  Great  Basin).  The  differences  in  part  reflected  that  the  diversity  and,  in  most 
cases,  the  complexity  of  vegetation  in  the  Chesapeake  Bay  Lowlands  was  greater  than  in  the 
Great  Basin.  Also,  our  methods  for  the  Great  Basin  were  consistent  with  the  methods  we  used 
during  previous  research  on  breeding  birds  in  the  region,  which  allowed  us  to  augment  the  data 
collected  for  this  project.  Although  field  methods  differed,  the  fundamental  information  captured 
in  both  regions  was  similar,  facilitating  inclusion  of  similar  covariates  in  models  and  comparison 
of  model  outputs  among  regions. 
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To  characterize  local  vegetation  composition  and  structure  at  the  bird- survey  points  in  the 
Chesapeake  Bay  Lowlands,  we  measured  four  radial  15-m  lines,  oriented  in  each  of  the  cardinal 
directions,  from  the  center  of  the  point  (Figure  6).  At  each  1  m,  we  recorded  the  identity  of  all 
plants  that  touched  each  10-cm  vertical  increment  of  a  1-m  tall  pole,  and  the  highest  point  on  the 
pole  at  which  each  plant  touched.  We  recorded  the  species  (if  possible)  or  the  functional  group  of 
all  saplings  <  1  cm  diameter  and  the  presence  of  shrubs,  non-native  Japanese  stiltgrass 
0 Microstegium  vimineum ),  other  grasses,  and  forbs.  We  also  recorded  whether  woody  debris  >  5 
cm  diameter,  leaf  litter,  water,  or  bare  ground  intersected  the  pole  on  the  ground.  We  used  a 
moosehorn  to  determine  whether  canopy  cover  at  each  1  m  was  nonzero. 

We  recorded  each  tree  (1-10  cm  diameter)  within  each  quadrant  (northeast,  southeast,  southwest, 
northwest)  of  a  circle  with  7.5-m  radius  extending  from  the  center  of  the  point.  Within  a  circle 
with  15-m  radius  extending  from  the  center  of  the  point,  we  also  recorded  the  genus,  diameter  at 
breast  height  (dbh),  and  stratum  (canopy  or  subcanopy)  of  all  live  trees  and  snags  >  10  cm  dbh. 
Within  the  latter  circle,  we  measured  the  dbh  (if  upright)  and  length  of  all  downed  wood  >10  cm 
dbh.  Limbs  of  branched  trees  were  measured  separately,  and  only  the  portion  of  the  wood  within 
the  circle  was  measured.  To  measure  understory  cover,  we  counted  the  number  of  decimeter 
segments  of  a  3-m  tall  pole  that  were  >75%  visible  when  the  pole  was  placed  15  m  from  the 
center  of  the  point  in  each  of  the  cardinal  directions. 

Vegetation  data  associated  with  birds  in  this  region  were  included  in  analyses  of  hierarchical 
occupancy  and  responses  of  songbird  occupancy  to  environmental  change. 


Figure  6.  Design  of  vegetation  sampling  for  birds  in  the  Chesapeake  Bay  Lowlands. 
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To  characterize  local  vegetation  composition  and  structure  at  the  bird- survey  points  in  the  central 
and  western  Great  Basin,  we  measured  three  radial  30-m  lines  from  the  center  of  the  point 
(Figure  7).  Lines  were  separated  from  each  other  by  120  degrees.  The  distal  end  of  each  line 
became  the  center  of  a  circular  vegetation  sampling  unit  with  1 1 .3-m  radius  (0.04  ha).  Within 
each  unit,  we  recorded  identities  and  sizes  of  all  live  trees  (either  dbh  or  basal  diameter, 
depending  on  plant  morphology).  We  recorded  the  identities  and  sizes  of  standing  dead  trees.  We 
used  a  concave  spherical  densiometer  to  estimate  proportion  of  canopy  cover.  To  estimate 
frequency  of  shrubs  and  ground  vegetation,  we  used  an  ocular  tube  with  measurements  taken  at  a 
45-degree  angle  downward  from  the  line  of  sight  (Noon  1981).  We  recorded  occurrence  of 
dominant  tree,  shrub,  and  herbaceous  taxa  (approximately  20-30).  We  collected  21  densiometer 
and  ocular  tube  readings  at  each  plot:  one  each  at  8  m,  16  m,  and  24  m  along  the  30-m  line  from 
the  center  of  the  plot  to  the  perimeter  of  each  circle,  and  one  while  facing  in  each  of  the  four 
cardinal  directions  from  the  center  of  each  circle.  Data  on  vegetation  and  other  covariates 
associated  with  birds  in  the  central  and  western  Great  Basin  were  included  in  analyses  of 
hierarchical  species  richness,  estimation  of  species  richness  of  multiple  taxonomic  groups  on  the 
basis  of  a  single  group,  and  partitioning  drivers  of  beta  diversity. 


Figure  7.  Design  of  vegetation  sampling  for  birds  in  the  central  and  western  Great  Basin. 


As  covariates  in  analyses  of  butterflies  in  the  Chesapeake  Bay  Lowlands,  we  measured  the 
length  (km)  of  all  edges  between  forest  and  agriculture,  ruderal,  or  herbaceous-developed  land 
cover;  structural  heterogeneity  of  the  understory  from  0-3  m  above  ground  (the  approximate 
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height  of  the  white-tailed  deer  browse  line  [Allombert  et  al.  2005,  Bressette  et  al.  2012]);  the 
proportion  of  the  basal  area  of  trees  (>  10  cm  dbh)  that  was  deciduous;  the  number  of  deciduous 
stems  (single-  or  multiple- stemmed  trees  or  shrubs;  1  to  <10  cm  dbh)  below  the  canopy;  and  the 
categorical  abundance  of  nectar.  We  included  categorical  abundance  of  nectar  as  a  detection 
covariate.  To  derive  edge  length,  we  first  obtained  data  on  land  cover  at  30-m  resolution  (2013 
Existing  Vegetation  Type  data;  www.landfire.gov).  Next,  we  delineated  edges  between  land- 
cover  types  in  Geospatial  Modeling  Environment  (www.spatialecology.com/gme/index.htm). 

We  then  derived  edge  length  in  ArcGIS  10.1  (ESRI,  Redlands,  California)  as  the  mean  of  the  30- 
m  cells  within  a  90-m  buffer  on  either  side  of  the  transect.  We  used  light  detection  and  ranging 
(lidar)  data  that  were  captured  from  22  April  -  10  May  2010  and  21-31  March  2013  to  estimate 
structural  heterogeneity  of  the  understory  on  the  basis  of  density  of  returns  at  10-m  resolution, 
averaged  among  the  10-m  cells  within  the  90-m  buffer.  We  measured  the  proportion  of  basal 
area  of  trees  that  was  deciduous  within  three  circular  plots  (15-m  radius)  that  were  randomly 
placed  within  the  90-m  buffer.  We  counted  the  number  of  deciduous  stems  within  three  circular 
plots  (7.5-m  radius),  each  of  which  was  embedded  within  one  of  the  15-m  plots.  Covariates 
associated  with  butterflies  in  the  Chesapeake  Bay  Lowlands  were  included  in  analyses  of 
occupancy  of  butterflies. 

As  covariates  in  analyses  of  butterflies  in  the  Great  Basin,  we  included  elevation,  the  square  of 
elevation,  terrain  roughness,  precipitation  in  the  water  year  (1  October  -  30  September)  of 
sampling,  and  categorical  abundance  of  nectar  and  mud  as  covariates  of  occupancy  in  the  Great 
Basin.  All  reasonably  might  be  expected  to  affect  habitat  quality  for  many  butterfly  species  (e.g., 
Fleishman  et  al.  2001a,  b).  We  included  categorical  abundance  of  nectar  and  mud  as  detection 
covariates.  We  derived  mean  elevation  of  the  transect  from  a  10-m  digital  elevation  model 
(www.ned.usgs.gov),  assuming  that  the  sampled  area  included  25  m  on  either  side  of  the 
transect.  We  used  a  digital  elevation  model  to  derive  terrain  ruggedness  (Riley  et  al.  1999)  within 
30-m  circular  neighborhoods  and  then  averaged  terrain  ruggedness  for  the  transect.  We  derived 
precipitation  at  4- km  resolution  from  the  Parameter-elevation  Relationships  on  Independent 
Slopes  Model  (PRISM).  Covariates  associated  with  butterflies  in  the  central  and  western  Great 
Basin  were  included  in  analyses  of  hierarchical  species  richness,  periodicity  of  sampling  (trade¬ 
offs  between  single-day  and  mutiple-days  sampling),  estimation  of  species  richness  of  multiple 
taxonomic  groups  on  the  basis  of  a  single  group,  occupancy  of  butterflies  in  diverse 
biogeographic  regions,  and  partitioning  drivers  of  beta  diversity. 

Response  of  songbird  occupancy  to  environmental  change.  We  measured  deer  use  and 
densities  of  forest  songbirds  in  two  areas  of  Virginia.  Our  primary  study  region,  coastal  Virginia, 
included  the  Virginia  and  Middle  Peninsula  on  the  coastal  plain  of  southeastern  Virginia  from 
Newport  News  to  West  Point.  In  this  area,  we  established  92  points  in  that  we  visited  three  times 
annually  from  2010-2013.  We  randomly  placed  these  points  in  wooded  areas  that  we  stratified 
on  the  basis  of  land  ownership  (private,  city,  county,  state,  and  federal).  Deer  appeared  to  be 
highly  abundant  in  coastal  Virginia.  Our  second  study  region,  inland  Virginia,  was 
approximately  200  km  west  of  the  coastal  Virginia  region,  and  encompassed  the  Shenandoah 
River  Valley  from  Stuarts  Draft  to  Harrisonburg.  In  this  area,  we  established  99  points  that  we 
surveyed  three  times  in  2012.  We  selected  sites  opportunistically  on  the  basis  of  access,  with  the 
goal  of  even  spatial  representation.  Deer  appeared  to  be  considerably  less  abundant  than  in 
coastal  Virginia. 
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We  used  distance  sampling  to  quantify  deer  pellets  and  forest  songbirds.  We  surveyed  birds  with 
8-min  variable-distance  point  counts.  Following  each  avian  survey,  we  counted  deer  fecal  pellets 
along  two  60-m  transects  randomly  placed  within  150  m  of  the  point  (coastal  Virginia)  or 
centered  on  the  point  (inland  Virginia).  We  assumed  that  pellet  density  reflected  the  relative 
intensity  of  deer  use.  Following  McShea  and  Rappole  (2000),  we  used  Partners  in  Flight  Species 
Assessment  Database  scores  for  bird  conservation  region  27  (Southeastern  Coastal  Plain)  as  a 
measure  of  conservation  priority  (Partners  in  Flight  Science  Committee  2012).  We  used  the 
regional  combined  score  for  the  breeding  season,  which  incorporates  the  species’  breeding 
distribution,  population  size,  regional  population  trend,  relative  breeding  density,  and  regional 
threats  to  reproduction.  Higher  scores  indicate  higher  conservation  priority. 


Analytical  methods 

Hierarchical  estimates  of  species  richness 

These  data  are  hierarchically  structured:  points  at  which  birds  were  surveyed,  or  transects  at 
which  butterflies  were  surveyed,  are  nested  within  canyons  (which  in  turn  are  nested  within 
mountain  ranges)  in  both  the  central  and  western  Great  Basin.  Our  response  variable  was 
cumulative  species  richness  across  all  years  in  which  the  location  was  surveyed.  We  modeled 
two  levels  of  the  data:  points  or  transects  and  canyons.  We  treated  canyon  as  a  random  effect, 
and  used  number  of  survey  years  as  a  control  variable.  We  also  estimated  values  of  both  point- 
level  covariates  and  canyon-level  covariates,  reflecting  our  hypothesis  that  species  richness  may 
respond  to  phenomena  that  are  relevant  at  different  spatial  resolutions.  We  ensured  that  highly 
collinear  covariates  were  not  included  in  the  same  analyses. 

In  the  formulae  below,  i  indexes  point  and  j  indexes  canyon.  x<  references  point-level  covariate  i, 
whereas  p,  references  canyon-level  covariate  j.  In  general,  when  we  treated  canyon  as  an 
indicator  variable,  we  could  apply  three  types  of  linear  regression  models: 

•  Constant  slopes  and  intercept:  y  =  a  +  2>, 

•  Constant  slopes,  varying  intercepts:  y  =  or + 

•  Varying  slopes  and  intercepts:  y  =  aj{i]  +  E/v. 

We  first  explored  an  unconditional  means  model,  which  does  not  include  covariates.  In  this 
model,  a  different  intercept  (random  effect)  is  fit  for  each  canyon.  This  hierarchical  model  is  fit 
to  estimate  components  of  variance  across  the  two  levels  (canyons  and  points).  The  model  is 

T2  2 

p  =  — - - ,  where  x  is  variance  among  the  random  intercepts  (i.e.,  mean  among-canyon 

r  +  co¬ 
variance),  and  cr2  is  residual  variance  (i.e.,  mean  variability  within  canyons).  We  used  the 
estimates  of  r2  and  cr2  to  compute  the  intraclass  correlation  coefficient,  p  ,  which  ranges  from 
0  (no  variance  among  canyons)  to  1  (no  within-canyon  variance). 


30 


Second,  we  fit  two  models  that  included  point-level  covariates  only.  The  first  of  these,  a  standard 
multiple-regression  model  with  fixed  effects,  had  constant  slopes  and  intercept  (a).  We  pooled 
data  across  canyons.  The  model  is  5,.  =  a  +  ^  j3ixi  +  st .  The  second  of  these  models  has  constant 

i 

slope  but  random  or  varying  intercept,  which  reflects  the  possibility  that  relations  between 
cumulative  species  richness  and  point-level  covariates  vary  among  canyons.  We  treated  the 
intercept  as  a  random  effect,  thereby  assuming  that  the  intercept  terms  had  a  normal  distribution. 
This  model  is 

Si=am+H^ixi+£i 

i 

a,  N^a^a) 


Third,  we  fit  two  models  that  included  canyon-level  covariates.  The  first  of  these  had  constant 
slopes  and  a  random  intercept.  The  intercept  could  be  a  function  of  the  canyon-level  covariates, 
Uj.  This  model  is 

Si=am+Y,Pixi+ei 

i 

a  -  a  ■  buj  +  rjj 

The  second  of  the  models  that  included  canyon-level  covariates  had  random  slopes  and  a  random 
intercept.  Both  the  slopes  and  the  intercept  could  vary  among  canyons,  as  a  function  of  the 
canyon-level  covariates.  This  model  is 

=  aj\i ]  +  '^jfij[i]Xi  +£i 

i 

aj=ao+b0uj+rjjl 

Pj=al+bluj+T)jl 


Periodicity  of  sampling:  duration  of  point  counts  of  birds 

To  estimate  detection  probabilities  (the  probability  of  detecting  a  species  at  a  site  if  it  is  present) 
(p)  and  occupancy  (the  probability  that  a  given  location  is  occupied  by  a  species)  (if/)  of  breeding 
birds  in  the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin,  we  used  a 
hierarchical,  single-season  (single-year)  occupancy  model  (MacKenzie  et  al.  2006).  We  modified 
the  general  form  of  the  model  to  allow  detection  probability  to  vary  among  observers 
(Diefenbach  et  al.  2003,  Alldredge  et  al.  2007). 

We  modeled  single-season  occupancy  (MacKenzie  et  al.  2002)  in  2012  and  2013  for  species 
detected  during  5-min  counts  (i.e.,  8-min  counts  truncated  at  5  min)  and  species  detected  during 
8-min  counts.  We  estimated  occupancy  for  all  species  with  a  detection  probability  >  0.3  in  both 
5-min  and  8-min  counts  in  a  given  year  and  for  which  the  confidence  intervals  around  the 
detection  probability  did  not  range  from  zero  to  one.  We  tested  whether  the  occupancy  estimates 
for  each  species  were  significantly  different  (the  95%  confidence  interval  centered  on  the 
difference  between  the  5-min  and  8-min  occupancy  estimate  did  not  include  zero  [Schenker  and 
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Gentleman  2001])  when  based  on  5-min  versus  8-min  counts.  We  report  the  percentage  (mean  ± 
SE)  of  species  for  which  occupancy  models  converged  and  produced  estimates  (i.e.,  95% 
confidence  interval  did  not  range  from  zero  to  one  and  other  diagnostics  did  not  indicate 
estimation  problems),  and  for  which  detection  probabilities  based  on  both  count  durations  were  > 
0.3. 

To  examine  whether  inferences  about  occupancy  over  two  years  were  consistent  over  longer 
periods  of  time,  we  also  fit  single-season  occupancy  models  to  data  from  2009-2013  on 
American  Robins  ( Turclus  migratorius),  MacGillivray’s  Warblers  ( Geothlypis  tolmiei),  and 
Vesper  Sparrows  ( Pooecetes  gramineus )  in  the  central  Great  Basin.  These  species  collectively 
span  gradients  of  land-cover  associations,  local  rarity,  and  ease  of  identification.  We  examined 
difference  in  the  precision  of  occupancy  estimates  based  on  5-min  and  8-min  counts.  We  defined 
precision  as  the  coefficient  of  variation  (CV).  We  conducted  all  analyses  in  package  Unmarked 
(Fiske  et  al.  2015)  in  R  (R  Core  Team  2014). 


Periodicity  of  sampling:  trade-offs  between  single-day  and  multiple-days  sampling 

We  fit  separate  models  of  occupancy  and  detection  (MacKenzie  et  al.  2002)  to  the  data  for  each 
species  (birds  in  the  Chesapeake  Bay  Lowlands,  anurans  in  the  Chesapeake  Bay  Lowlands,  and 
butterflies  in  the  central  and  western  Great  Basin),  sampling  design,  and  year.  For  all  species  we 
fit  models  that  represented  alternate  hypotheses  about  the  variation  in  detection  probability  (p ) 
among  sites  and  surveys,  but  did  not  use  covariates  to  model  variation  in  occupancy  (i//)  among 
sites.  We  fit  a  model  with  constant p  (i//  [.]  p  [.],  where  indicates  no  variation  in  the 
parameter),  and  a  model  with  survey-specific  pj  (i//  [.]  p  [t],  where  j  indexes  survey  and  is  treated 
as  a  categorical  variable  with  three  levels). 

For  birds,  we  also  fit  models  that  included  one  of  four  survey- specific  covariates:  time  (number 
of  minutes)  since  sunrise,  Julian  date,  observer  identity,  and  detection  history.  Detection  history 
refers  to  the  series  of  species-specific  detections  (represented  as  Is)  and  non-detections 
(represented  as  0s).  For  example,  the  detection  history  of  a  species  that  was  detected  on  the  first 
and  third  visits  but  not  detected  on  the  second  visit  would  be  101.  The  lag  between  sunrise  and 
initiation  of  a  survey  could  affect  p  because  song  frequency  decreases  as  time  since  sunrise 
increases  (Skirvin  1981).  We  included  Julian  date  to  capture  variation  in  p  throughout  the 
breeding  season  (Best  1981,  Skirvin  1981).  We  included  observer  identity  because  the  ability  to 
detect  a  given  species  varies  among  observers  (Bart  1985,  Diefenbach  et  al.  2003,  Alldredge  et 
al.  2007,  Simons  et  al.  2007,  Miller  et  al.  2011).  Inclusion  of  detection  history  allowed  estimates 
of  p  to  differ  between  sites  at  which  the  species  previously  had  been  detected  or  not  detected 
(MacKenzie  et  al.  2006).  Higher  estimates  of  p  at  the  former  suggest  dependence  among 
observations. 

For  anurans,  we  also  fit  models  that  included  one  of  four  survey-specific  covariates:  Julian  date, 
time  since  sunset,  temperature  (°C),  and  log  temperature  (a  pseudo-threshold  [Scherer  et  al. 
2012],  log  [°C]).  We  included  a  fifth  survey- specific  covariate,  number  of  days  since  above- 
average  rainfall,  for  species  that  breed  during  the  middle  and  late  breeding  season.  We  included 
Julian  date  to  account  for  temporal  variation  in  vocalization.  We  included  time  since  sunset 
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because  vocalization  frequency  may  vary  throughout  the  night.  We  included  temperature  because 
vocalization  frequency  may  decrease  as  temperature  decreases.  We  included  days  since  above- 
average  rainfall  because  species  that  breeding  during  the  middle  and  late  breeding  season 
vocalized  more  frequently  after  a  major  rain  event.  We  did  not  include  observer  identity  as  a 
covariate  because  there  was  some  consistency  in  observers  among  years.  We  did  not  include 
detection  history  as  a  covariate  because,  unlike  birds,  anurans  are  stationary  during  the  breeding 
season. 

For  butterflies,  we  fit  models  with  constant  and  survey- specific  p.  We  also  fit  models  that 
included  one  of  two  survey-specific  covariates,  relative  abundance  of  nectar  and  relative 
abundance  of  mud. 

We  derived  model-averaged  estimates  of  if/  (Burnham  and  Anderson  1998)  and  unconditional 
95%  confidence  intervals.  We  considered  differences  in  p  or  p *  (the  probability  that  a  given 
species  was  detected  at  least  once  during  three  surveys)  between  single-day  and  multiple-day  to 
be  statistically  significant  if  the  95%  confidence  interval  centered  on  the  difference  between  the 
single-  and  multiple-day  p  or  p*  estimate  did  not  include  zero  (Schenker  and  Gentleman  2001). 

p*  identifies  potential  reductions  in  precision  of  if/  estimates  between  sampling  designs  because 
naive  ip  is  divided  by  p*  to  estimate  if/  (MacKenzie  et  al.  2006).  As  p *  decreases,  if/  increases.  To 
calculate  p*,  we  first  compared  model- averaged  and  site-  and  survey-specific  p.  When  the  mean 
and  median  of  site-  and  survey-specific  estimates  of  p  were  similar,  we  used  estimates  of  p  from 
the  model  with  constant  p  to  estimate  p*. 

We  examined  whether  the  output  of  each  model  included  evidence  of  lack  of  convergence  or 
estimation  problems  (e.g.,  extremely  large  standard  errors).  If  we  observed  either  of  the  latter  in 
at  least  one  model  for  a  given  species  in  either  survey  year,  we  removed  the  species  from  further 
analyses  for  that  year.  In  other  words,  we  restricted  our  inferences  to  those  species  for  which  all 
models  converged  in  both  designs  in  one  year. 

We  report  the  percentage  (mean  ±  SE)  of  species  for  which  occupancy  models  converged  and 
produced  reasonable  estimates  (i.e.,  95%  confidence  interval  did  not  range  from  zero  to  one  and 
other  diagnostics  did  not  indicate  estimation  problems).  We  fit  all  models  with  R  Mark  (Laake 
2013)  in  the  R  statistical  software  package  (R  version  2.15.2,  http:/www.r-project.org/,  accessed 
15  Dec  2013). 

Simulations.  To  evaluate  potential  biases  inherent  in  single-day  or  multiple-days  sampling 
designs,  we  simulated  data  in  which  detection  probability  and  occupancy  were  constant  but  the 
species’  availability  varied.  We  based  known,  taxonomic  group-specific  detection  probabilities 
and  occupancy  on  the  average  estimate  across  all  species  for  which  models  converged.  The 
availability  process  is  random  in  the  multiple-days  design.  However,  the  availability  process  is 
Markovian  in  the  single-day  design;  if  a  given  species  is  available  during  the  first  survey,  it  is 
more  likely  to  be  available  during  the  second  and  third  survey.  As  a  result,  heterogeneous 
detection  histories  (e.g.,  010,  100,  101)  are  more  common  in  data  from  the  multiple-days  design 
and  homogeneous  detection  histories  (e.g.,  000,  111)  are  more  common  in  data  from  the  single¬ 
day  design. 
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For  each  design,  we  varied  availability  in  0.1  increments  from  0.1  to  1.0.  To  estimate  v| /  and  p* 
for  each  of  the  ten  availability  values,  we  first  created  detection  histories  for  130  sites  and 
repeated  this  process  100  times.  Presence  or  absence  was  random  for  each  survey  in  the  multiple- 
days  design  and  the  first  survey  in  the  single-day  design.  The  second  and  third  surveys  in  the 
single-day  design  had  the  same  state  as  the  first  survey.  For  each  of  100  detection  histories  we 
then  estimated  \\i  and  p*  and  derived  a  final  mean  and  standard  deviation  for  each  availability 
value  on  the  basis  of  the  100  runs.  We  ran  all  models  with  R  Mark  (Laake  2013)  in  the  R 
statistical  software  package  (R  version  2.15.2,  http:/www.r-project.org/,  accessed  15  Dec  2013). 


Estimation  of  species  richness  of  multiple  taxonomic  groups  on  the  basis  of  a  single  group 

Matching  of  survey  data.  Several  of  our  analyses  required  comparisons  of  point-based  surveys 
of  birds  with  transect-based  surveys  of  butterflies.  To  generate  spatially  comparable  survey  units, 
we  buffered  butterfly  transects  by  100  m  (the  radius  of  the  bird  points)  on  either  side  and 
compiled  records  of  all  bird  species  recorded  at  points  within  the  buffered  butterfly  transect.  We 
excluded  all  bird  points  that  did  not  fall  within  a  buffered  butterfly  transect  and  all  butterfly 
transects  that  included  no  bird  points.  Differences  in  the  years  in  which  birds  and  butterflies  were 
surveyed  led  to  differences  among  models  with  respect  to  the  years  of  data  that  were  used  for 
indicator  selection  and  evaluation  (see  below). 

Environmental  data.  We  evaluated  whether  1 1  environmental  variables — canyon  area, 
elevation,  terrain  ruggedness,  minimum  temperature,  maximum  temperature,  precipitation, 
prevalence  of  canopy,  prevalence  of  shrubs,  number  of  large  trees,  riparian  fragmentation,  and 
riparian  cover — explained  or  predicted  the  species  richness  of  birds.  We  selected  these  variables 
on  the  basis  of  their  relevance  to  the  Great  Basin  avifauna  (Dobkin  and  Wilcox  1986,  Ehrlich  et 
al.  1988)  We  evaluated  seven  environmental  variables  relative  to  the  species  richness  of 
butterflies  that  one  reasonably  might  expect  to  be  relevant  to  that  taxonomic  group  in  the  Great 
Basin  (Boggs  and  Ross  1993,  Fleishman  et  al.  2003):  sampled  area,  mean  elevation,  terrain 
ruggedness,  maximum  temperature,  precipitation,  abundance  of  nectar,  and  abundance  of  mud. 

To  estimate  canyon  area  (ha),  we  connected  the  centroids  of  the  points  along  the  minor  road  or 
trail  in  each  canyon  (on  or  near  which  we  sampled)  and  buffered  100  m  on  either  side  of  that 
line.  We  calculated  elevation  (m),  which  we  derived  from  a  10-m  digital  elevation  model 
(DEM),  as  the  mean  elevation  within  100  m  of  the  center  of  the  point  (birds)  or  within  the 
transect  buffered  25  m  on  either  side  (butterflies).  We  calculated  terrain  ruggedness  (Riley  et  al., 
1999),  which  we  also  derived  from  a  10-m  DEM,  as  the  mean  within  the  100-m  point  (birds)  or 
within  30-m  circular  neighborhoods  over  the  transect  (butterflies).  We  derived  the  three  weather 
variables  at  4-km  resolution  from  the  Parameter-elevation  Relationships  on  Independent  Slopes 
Model  (PRISM)  for  the  water  year  (October-September)  in  which  a  given  survey  was  conducted. 
We  calculated  maximum  and  minimum  temperature  as  the  mean  of  mean  maximum  (or 
minimum)  monthly  temperature  (°C),  and  precipitation  as  the  sum  of  the  monthly  precipitation 
totals.  We  measured  prevalence  of  canopy  and  shrubs  as  the  proportion  of  the  21  sampling 
locations  in  a  given  point  at  which  canopy  or  shrubs  was  present.  Trees  were  classified  as  large  if 
their  diameter  at  breast  height  (dbh)  or  basal  diameter  was  >16  cm.  We  calculated  riparian 
fragmentation  at  the  extent  of  the  canyon  (30-m  resolution)  as  the  mean  distance  from  each 
riparian  cell  to  the  nearest  riparian  cell  within  the  sampled  area  of  the  canyon.  Proportion  of 
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riparian  cover  within  each  point  was  derived  from  measures  of  the  normalized  difference 
vegetation  index  (NDVI)  (generally  >0.25)  in  National  Agriculture  Imagery  Program  (NAIP) 
aerial  images  (https://lta.cr.usgs.gov/NAIP). 

We  used  year-specific  values  for  variables  that  were  highly  dynamic  and  for  which  we  had 
multiple  years  of  data  (temperature,  precipitation,  and  abundance  of  nectar  and  mud).  The 
topographic  variables  we  measured  were  static  over  the  time  period  of  our  work.  Values  of  our 
vegetation  variables  unlikely  to  change  substantially  over  less  than  a  decade  in  the  absence  of 
major  disturbances. 

Indicator  selection.  We  used  Bayesian  model  selection  to  identify  sets  of  five  indicator  species 
that  explained  the  greatest  proportion  of  variation  in  species  richness  ( indicator  selection).  Our 
primary  focus  was  the  development  of  a  transferable  method  for  identifying  indicator  species;  we 
did  not  necessarily  expect  that  the  same  set  of  species  would  be  associated  with  species  richness 
in  all  ecological  or  management  circumstances.  We  fitted  a  log-linear  Poisson  regression  model, 
but  estimated  model  coefficients  with  a  reversible -jump  Markov  chain  Monte  Carlo  sampler, 
which  allows  variables  to  be  removed  from  or  added  to  the  model  during  model  fitting  (Lunn  et 
al.  2008).  This  method  estimates  the  posterior  probability  of  inclusion  of  each  candidate 
indicator  species.  We  identified  the  five  species  with  the  highest  posterior  probabilities  of 
inclusion. 

The  response  variable,  species  richness  of  birds  or  species  richness  of  butterflies,  included  all 
species  detected  in  the  specified  subregion  or  time  period.  We  restricted  candidate  indicator 
species  to  those  detected  in  both  subregions  or  both  time  periods  because  only  these  can  serve  as 
transferable  indicators. 

Model  fitting.  We  fitted  random  forest  models  with  either  the  five  indicator  species,  all  of  the 
environmental  variables,  or  the  five  indicator  species  and  all  of  the  environmental  variables  (full 
models )  (Figure  8).  We  fitted  models  with  indicator  species  only  and  with  both  indicator  species 
and  environmental  variables  to  account  for  possible  confounding  between  indicator  species  and 
environmental  variables.  We  fitted  each  model  variant  ten  times  (hereafter  referenced  as  ten 
iterations).  We  report  results  as  averages  across  the  ten  iterations. 
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Model  fitting 


Model  evaluation 


Figure  8.  Main  steps  in  the  model  fitting  and  model  evaluation  process. 


The  fitted  models  included  a  different  value  of  the  response  variable  for  each  year.  For  example, 
if  we  built  a  model  with  data  collected  from  2012-2014,  the  model  included  three  values  of 
species  richness  for  each  location.  We  fitted  models  with  and  without  spatial  and  temporal 
random  effects  to  assess  whether  correlations  among  locations  or  years  affected  model  fit.  The 
inclusion  of  random  effects  substantially  reduced  the  accuracy  of  external  tests  of  the  fitted 
models.  Therefore,  we  report  results  for  models  fitted  without  these  terms. 

Model  evaluation.  We  used  the  three  types  of  fitted  models  to  estimate  species  richness  in  new 
locations  or  at  future  times  ( model  evaluation,  i.e.,  external  evaluation  or  prediction  [Elith  and 
Leathwick  2009];  in  the  case  of  the  indicator  species-only  models,  we  reference  indicator 
evaluation)  (Figure  8). 

We  conducted  separate  evaluations  of  the  spatial  and  temporal  transferability  of  our  fitted 
models  because  we  wished  to  avoid  confounding  the  two  types  of  transferability.  For  evaluations 
of  spatial  transferability,  we  conducted  indicator  selection  in  one  zoogeographic  subregion  (i.e., 
the  central  or  western  Great  Basin),  and  indicator  evaluation  in  the  other  subregion.  For 
evaluations  of  temporal  transferability,  we  conducted  indicator  selection  with  data  from  a  subset 
of  the  sampled  years  in  one  subregion,  and  indicator  evaluation  with  data  from  a  later  subset  of 
the  sampled  years  in  the  same  subregion.  For  each  model  variant,  we  used  each  of  the  ten  fitted 
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random  forest  models  to  predict  species  richness  in  the  other  subregion  or  for  the  later  time 
period.  We  assessed  model  accuracy  with  r2  values  (based  on  Pearson’s  r)  between  predicted  and 
observed  species  richness. 

We  present  and  discuss  the  number  and  identity  of  indicator  species  and  environmental  variables 
included  in  the  models  in  which  >0.25  of  the  variation  in  the  response  variable  was  predicted. 

We  used  this  threshold,  which  is  arbitrary  but  has  been  recommended  (Mpller  and  Jennions 
2002),  to  simplify  reporting  of  the  results.  Results  at  other  thresholds  (e.g.,  0,  0.50)  were  not 
substantially  different. 


Theory  and  practical  application  of  occupancy  models:  hierarchical  models  of  occupancy 

We  used  a  dynamic  occupancy  model  to  evaluate  associations  of  covariates  with  occupancy, 
persistence,  colonization,  and  detection  of  birds  in  the  Chesapeake  Bay  Lowlands  and  central 
and  western  Great  Basin  while  accounting  for  the  possibility  of  false  absences.  We  considered 
the  true  presence  or  absence  of  each  of  J  species,  each  of  which  was  a  member  of  one  of  G 
guilds,  at  each  of  K  points,  in  each  of  T  time  steps  to  be  a  partly  observed  state  variable  zj5k,t, 
which  is  equal  to  one  if  species  j  is  present  at  point  k  in  time  step  t,  and  equal  to  zero  otherwise. 
Because  we  used  a  model  specification  that  marginalized  over  the  latter  state  variable,  we 
directly  modeled  vpj,k,t,  the  probability  of  occurrence. 

Observation  model 

We  represented  the  detection  data  from  the  Chesapeake  Bay  Lowlands  as  yj,k,t,  which  indicates 
whether  species  j  was  detected  at  point  k  during  year  t.  Three  surveys  were  conducted  at  each 
point  in  each  year,  so  yj,k,t  was  a  vector  of  length  3  in  which  each  element  had  a  value  of  0  (not 
detected)  or  1  (detected).  Following  the  likelihood  formulation  of  Royle  and  Dorazio  (2008),  and 
marginalizing  the  latent  occupancy  state  variable,  our  likelihood  function  was 
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[Vj,k,t  I  n  Bernoulli^, fc(f>i 
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where  vpj,k,t  is  the  probability  of  species  j  occurring  at  point  k  in  year  t,  pj,k,i  is  the  probability  of 
detection  on  survey  i,  '  ^'=  Vj  k  '  is  an  indicator  function  that  takes  on  the  value  1  only  if 

_  q 

*"  V  *  '  and  is  zero  otherwise,  and  square  brackets  represent  the  probability  function 
(i.e.,  a  probability  mass  function  Pr(Y  =  y|0)  for  discrete  variables  or  a  probability  density 
function  p(y|0)  for  continuous  variables,  where  0  represents  parameters  of  the  distribution  of  the 
random  variable  Y. 


We  represented  the  detection  data  from  the  central  and  western  Great  Basin  as  yij,  which  is  the 
number  of  times  species  j  was  detected  during  year  i,  bounded  between  0  and  ni,  the  number  of 
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surveys  in  year  i.  Following  the  likelihood  formulation  of  Royle  and  Dorazio  (2008),  and 
marginalizing  the  latent  occupancy  state  variable,  our  likelihood  function  for  surveys  i  =  1, ...,  n 
and  species  j  =  1, ...,  J  was 


[yi,j  | '|/j ,kj ,t i ,pk i]=v|/j . kj ,t i B i nom iaUyi.j  lpki;ni)+I(yi j=0)( l~f j,ki,ti), 

where  xpj ,ki,ti  is  the  probability  of  species  j  occurring  at  point  ki  in  year  ti,  pj  is  the  probability  of 
detection,  I(yj  j  _q)  is  an  indicator  function  that  takes  on  the  value  1  only  if  yj  .j  =  0  and  is  zero 

otherwise,  and  square  brackets  represent  the  probability  function  (i.e.,  a  probability  mass 
function  Pr(Y  =  y|9)  for  discrete  variables  or  a  probability  density  function  p(y|9)  for  continuous 
variables,  where  9  represents  parameters  of  the  distribution  of  the  random  variable  Y  . 

Process  model.  The  process  model  has  three  components:  initial  occupancy,  colonization,  and 
extinction.  Colonization  and  extinction  determine  occupancy  beyond  the  initial  time  step.  For  all 
of  these  components,  we  created  a  set  of  covariates  that  reasonably  could  be  expected  to  be 
associated  with  occupancy  of  each  guild  (allowing  relations  to  vary  among  guilds). 

For  initial  occupancy  in  the  Chesapeake  Bay  Lowlands,  we  included  the  following  point-level 
covariates:  number  of  snag  saplings;  cover  of  leaf  litter;  volume  of  downed  woody  debris; 
number  of  deciduous  trees  in  the  subcanopy;  number  of  deciduous  trees  in  the  canopy;  number 
of  deciduous  saplings;  number  of  coniferous  saplings;  mean  canopy  height  within  100  m  of  the 
point  center;  structural  heterogeneity  between  0  and  3  m  height  within  100  m  of  the  point  center; 
proportion  of  riparian  cover  within  500  m  of  the  point  center;  proportion  of  mesic  deciduous  or 
mixed  tree  cover  within  500  m  of  the  point  center;  length  of  forest  edge  within  500  m  of  the 
point  center;  proportion  of  low  density,  moderate  density,  and  high  density  development  within 
500  m  of  the  point  center;  and  proportion  of  conifer  cover  within  500  m  of  the  point  center.  We 
embedded  this  information  in  a  site-level  design  matrix  Xs  in  which  the  columns  contained 
centered  and  scaled  continuous  covariates  that  interacted  with  guild.  This  design  matrix  therefore 
represented  the  interactions  of  three  guilds  with  16  covariates,  yielding  48  columns  and  K  x  G 
rows,  where  each  row  represented  a  site  by  guild  combination.  Because  occupancy  varies  among 
species,  we  included  a  species-specific  adjustment  for  initial  occupancy,  aj.  Therefore,  our 
submodel  for  initial  occupancy  was 

logit(v|/j,k,l)  =  px|/l  +aj^  +Xs(§J’k)pv|/l 

where  pivpi  is  a  global  intercept,  aj^  is  a  species  specific  adjustment,  Pv|/l  is  a  parameter  vector 

of  length  K  x  G,  and  Xs^J,k^  is  a  row  vector  from  Xs  corresponding  to  point  k  and  the  guild  of 
species  j(gj). 

For  subsequent  time  steps  t  >  1,  we  modeled  detection  probability  with  an  autoregressive 
formulation  that  included  a  persistence  probability  from  time  step  t  -  1  to  t  (cpj,k,t-l)  and  a 
colonization  probability  (the  probability  that  a  point  unoccupied  in  t  -  1  is  occupied  in  time  t): 
Yj,k,t-l: 

Vj,k,t  =  Yj,k,t-lcpj,k,t-l  +  (1  -  v|/j,k,t-l)yj,k,t-l 
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To  explore  potential  environmental  associations  with  colonization  and  extinction,  we  also 
modeled  (p  and  y  as  functions  of  the  same  point-level  covariates  introduced  for  the  initial 
occupancy  submodel.  We  also  included  species-specific  adjustments  through  random  effects  to 
account  for  variation  in  colonization  and  extinction  probabilities  among  species.  The  submodels 
for  persistence  and  colonization  therefore  were 

logit(cpj,k,t)  =  p  cp  +  Xs^jk)p  cps  + 
logit(yj,k,t)  =  py  +  Xs*Jk)pys  + 

where  and  c/(^  are  species-specific  adjustments  for  persistence  and  colonization,  P  terms 
represent  parameter  vectors,  and  p  terms  represent  global  intercepts. 

We  allowed  for  point-level  variation  in  detection  probabilities  by  modeling  the  effect  of  log- 
transformed  wind  speed,  time  since  sunrise,  and  days  since  May  1,  centered  and  scaled  to 
comprise  three  columns  in  a  detection  design  matrix  Xp,  which  had  K  x  T  x  3  rows,  one  for  each 
point  by  year  by  survey  combination.  We  also  included  a  species -specific  adjustment  in  detection 
probability: 

logit(pj,k,t  )=p  +Xp^')Pp 

For  initial  occupancy  in  the  central  and  western  Great  Basin,  we  included  the  following 
covariates:  mountain  range  (Shoshone,  Toiyabe,  Toquima,  and  Monitor  in  the  central  Great 
Basin;  Wassuk,  Sweetwater,  and  Sierra  in  the  western  Great  Basin),  canopy  prevalence  (the 
number  of  presences  divided  by  the  number  of  samples,  generally  21  per  point),  elevation  (mean 
elevation  within  the  center  of  the  point,  derived  from  a  10-m  digital  elevation  model),  mean 
canopy  cover,  mean  coniferous  canopy  cover,  mean  riparian  canopy  cover,  rabbitbrush  incidence 
(number  of  presences  divided  by  number  of  samples),  riparian  fragmentation  (calculated  as  the 
mean  distance  from  each  riparian  cell  to  the  nearest  riparian  cell  within  the  sampled  area), 
sagebrush  prevalence  (number  of  presences  divided  by  number  of  samples),  and  terrain 
ruggedness.  We  embedded  the  covariate  information  in  a  point-level  design  matrix  Xs  in  which 
the  columns  contained  indicator  variables  for  categorical  covariates  (such  as  mountain  range), 
and  centered  and  scaled  continuous  covariates,  both  interacting  with  guild. 

The  full  design  matrix  represented  the  interactions  of  three  guilds  with  each  of  the  12  covariates 
(36  columns)  and  K  x  G  rows,  where  each  row  represents  a  point-by-guild  interaction.  Because 
occupancy  varies  among  species,  we  included  a  species-specific  adjustment  for  initial 
occupancy,  aj.  Our  submodel  for  initial  occupancy  was 


logit(v|/j,k,l)=pv|/i  +  aj  +Xs(gJ’k)pi|/l, 

where  pivpi  is  a  global  intercept,  aj  is  a  species  specific  adjustment,  f3\|/  ]  is  a  parameter  vector  of 
length  K  x  G,  and  Xs  (gj  ,k)  is  a  row  vector  from  Xs  corresponding  to  point  k  and  the  guild  of 
species  j(gj). 

For  subsequent  time  steps  t  >  1,  we  modeled  detection  probability  with  an  autoregressive 
formulation  that  included  a  persistence  probability  from  time  step  t  -  1  to  t:  (pj,k,t-F  and  a 
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colonization  probability  (the  probability  that  a  point  unoccupied  in  t  -  1  is  occupied  in  time  t): 
Yj,k,t— 1 ; 


Yj,k,t  =  x|/j,k,t-l(pj,k,t-l  +  (1  -  Yj,k,t-l)Yj,k,t-l. 

To  explore  spatial  and  temporal  associations  with  colonization  and  extinction,  we  also  modeled 
cp  and  y  as  functions  of  covariates. 

For  the  central  and  western  Great  Basin,  we  included  both  the  same  point-level  covariates 
introduced  for  the  initial  occupancy  submodel  and  a  new  set  of  spatially  and  temporally  dynamic 
covariates  at  the  year  level:  maximum  air  temperature,  minimum  air  temperature,  and  total 
precipitation,  all  from  March  through  May.  These  three  covariates  crossed  with  the  three  guilds 
generated  9  variables  that  comprised  the  columns  of  a  dynamics  design  matrix  Xd,  which  had  K 
x  G  x  (T  -  1)  rows  (among  T  time  steps,  T  -  1  colonization  and  persistence  events  are  possible). 
The  submodels  for  persistence  and  colonization  thus  were 

logit(cp  j,k,t)  =  pep  +  Xs^pcps  +  Xd^J  ’^p  ^ 
logit(y  j,k,t)  =  py  +  Xs^)p  ys  +  Xd^J  ’^p  yd 

where  Xd^^  extracts  the  row  from  Xd  corresponding  to  the  guild  of  species  j  (gj),  point  k,  and 
time  step  t,  P  terms  represent  parameter  vectors,  and  p  terms  represent  global  intercepts. 

We  allowed  for  point-level  variation  in  detection  probabilities  by  modeling  the  effect  of  canopy 
cover,  shrub  cover,  and  terrain  ruggedness,  centered  and  scaled  to  comprise  three  columns  in  a 
detection  design  matrix  Xp,  which  had  K  rows,  one  for  each  point: 

logit(pk)  =  pp  +  Xp^)pp 

Parameter  model.  To  complete  our  model  specifications  for  the  Chesapeake  Bay  Lowlands  and 
central  and  western  Great  Basin,  we  specified  prior  distributions  for  all  remaining  parameters  as 
follows.  All  of  the  global  intercepts  (ppsil,pphi,py)  received  Normal(p  =  0,a  =  1.5)  priors,  which 
led  to  nearly  uniform  prior  distributions  on  the  probability  scale.  The  species-specific 
adjustments  for  initial  occupancy  were  treated  as  normally  distributed  random  effects  with 
unknown  variance:  a  E  Normal(0,ca),  with  o<x  E  Normal(0,  1).  For  our  three  detection 
coefficients,  we  assigned  conservative  priors:  (3p  E  Normal(0,  1). 

Across  the  initial  occupancy,  persistence,  and  colonization  submodels,  we  had  well  over  100  P 
coefficients  to  estimate.  However,  we  expected  a  priori  that  the  majority  of  these  coefficients 
would  be  nearly  zero.  Therefore,  we  used  a  scale  mixture  prior,  adapted  from  the  horseshoe 
prior,  which  induces  sparseness  by  shrinking  estimates  toward  zero  while  allowing  for  some  of 
the  coefficients  to  be  far  from  zero.  Concatenating  [h|/] ,  Pcps,  P(pcp  Pys,  and  Py^  into  one  vector 

P  ,  we  assumed  the  following  scale  mixture  prior  for  parameters  1  =  1, ...,  L,  where  L  is  the  total 
number  of  coefficients  in  P  : 
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Pl*0  Normal(0,x^l) 

Xl  0  Student-t(0,  1,  r|  =  3) 
x  0  Normal(0,  1) 

This  prior  distribution  is  similar  to  the  original  horseshoe  prior,  but  uses  a  lighter  tailed  Student-t 
distribution  for  X  rather  than  a  Cauchy  distribution,  which  is  more  reasonable  under  the  logit  link 
function. 


Use  of  occupancy  as  a  surrogate  measure  of  species  richness:  estimation  of  the  occupancy 
of  butterflies  in  diverse  biogeographic  regions 

We  used  single-season  occupancy  models  to  analyze  the  data  collected  from  all  transects  and  on 
all  visits  during  2013  and  2014  for  13  species  of  butterflies  in  the  Chesapeake  Bay  Lowlands  and 
15  species  each  in  the  central  and  western  Great  Basin.  We  generally  restricted  our  analyses  to 
species  that  were  detected  in  >  30%  and  <  70%  of  the  transects  in  each  year  and  that  are  not 
migrants,  highly  vagile  (e.g.,  thousands  of  meters),  or  do  not  complete  their  entire  life  cycle  in 
the  ecosystem.  We  also  modeled  single-season  occupancy  of  the  15  species  in  the  central  Great 
Basin  in  1995. 

We  fit  two  parameterizations  (MacKenzie  et  al.  2002,  Kendall  et  al.  2013)  of  the  single-season 
occupancy  model.  Both  include  the  parameters  y/i,  the  probability  that  a  given  species  occupies 
transect  i,  and  prh  the  probability  that  the  species  is  detected  given  that  it  is  present  on  transect  i 
during  visit  j.  The  Kendall  et  al.  (2013)  model  also  allows  a  single  entry  and  exit  of  the  species 
from  each  transect  during  the  sampling  period;  the  probabilities  of  entry  and  exit  between  visits  j 
and  j  +  1  are  denoted  as  pj  and  dij,  respectively.  py  may  vary  during  the  period  in  which  the 
species  is  available. 

We  fit  models  to  the  occupancy  data  in  two  stages.  In  the  first  stage,  we  evaluated  submodels  of 
py,  pij ,  and  dij,  and  tested  the  assumption  of  closure.  For  both  parameterizations,  the  sets  of 
submodels  of  py  included  effects  of  categorical  abundance  of  nectar,  mud  (in  the  Great  Basin),  or 
both.  In  all  cases,  we  used  the  highest-abundance  class  of  nectar  or  mud  as  the  intercept. 

We  also  included  a  fixed  effect  of  visit  (i.e.,  we  estimated  py  for  each  survey)  in  the  submodels. 
Furthermore,  the  submodels  included  additive  effects  of  nectar  and  visit  and  of  mud  and  visit. 

For  the  models  in  which  the  closure  assumption  was  relaxed,  we  estimated  ft]  and  dij  as  linear 
and  quadratic  functions  of  visit.  We  specified  full  models  that  included  every  combination  of  the 
submodels  of  pij,  Py,  and  dij  with  an  intercept-only  submodel  of  occupancy.  We  used  Akaike’s 
Information  Criterion  adjusted  for  small  sample  sizes  (AICc)  and  Akaike  weights  (wm),  where  m 
indexes  models,  to  compare  submodels  of  pij,  pij,  and  dij  (Burnham  and  Anderson  2002).  We 
retained  all  models  from  the  first  stage  of  modeling  with  AICc  values  within  2  units  of  the 
highest-ranked  model  and  included  them  in  the  second  stage  of  modeling. 

In  the  second  stage  of  modeling,  we  evaluated  submodels  of  y/i.  The  set  of  submodels  included 
effects  of  each  covariate  and  several  ecosystem-specific  interactions.  We  combined  all 
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submodels  of  y/i  with  the  submodels  of  other  parameters  from  the  model  that  was  ranked  highest 
in  the  first  stage.  Because  many  models  that  included  covariates  generated  highly  imprecise 
occupancy  estimates,  we  report  y/i  from  the  intercept- only  submodel. 

We  standardized  and  centered  all  continuous  covariates.  We  calculated  Pearson  product-moment 
correlations  between  continuous  covariates.  We  did  not  include  two  continuous  covariates  in  a 
given  model  if  their  correlation  coefficient  was  >  0.60.  We  examined  box-and-whisker  plots, 
created  in  R  (R  Core  Team,  2013),  to  assess  correlations  between  continuous  and  categorical 
covariates.  We  used  the  plots  to  anticipate  potential  confounding  effects  of  multicollinearity. 
When  a  strongly  supported  model  included  both  a  continuous  and  a  categorical  covariate,  we 
examined  whether  the  magnitude  or  direction  of  regression  coefficients  in  the  latter  model  and  in 
models  that  included  each  of  those  covariates  alone  was  considerably  different.  We  included  a 
maximum  of  two  covariates  in  additive  models. 

We  characterized  the  strength  of  association  between  response  variables  and  covariates  on  the 
basis  of  the  AICc  values  of  the  models  in  which  they  were  included  and  the  degree  to  which 
estimates  of  the  95%  confidence  intervals  (CIs)  of  the  regression  coefficients  overlapped  zero.  If 
a  covariate  was  included  in  the  model  with  the  lowest  AICc,  or  in  a  model  with  an  AICc  value 
within  2  units  of  the  model  with  the  lowest  AICc,  we  considered  it  to  be  associated  with  /?</  or  y/i 
and  report  it  here.  We  considered  the  strength  of  association  of  a  covariate  with  py  or  y/i  to  be 
greater  if  its  CIs  did  not  overlap  zero  than  if  its  CIs  overlapped  zero.  We  report  associations  with 
continuous  covariates  as  regression  coefficients  and  associations  with  categorical  covariates  as 
effect  sizes. 

When  data  are  limited  or  probabilities  approach  0  or  1,  parameters  may  not  be  estimated 
correctly.  Evidence  of  incorrect  estimates  includes  noticeably  high  values  of  parameters  or  their 
standard  errors  and  estimates  of  standard  errors  that  are  near  zero.  We  examined  estimates  of 
model  parameters  and  used  the  diagnostics  in  Program  MARK  to  identify  potentially 
questionable  estimates.  We  also  considered  regression  coefficients  or  effect  sizes  with  absolute 
values  >10  to  be  questionable. 


Use  of  occupancy  as  a  surrogate  measure  of  species  richness:  environmental  associations 
with  multiple  states  of  anuran  occupancy 

For  each  species,  we  combined  the  data  from  the  wetland  surveys  in  201 1  and  2012  and  fit  both 
a  single-season,  single-state  occupancy  model  (MacKenzie  et  al.  2002)  and  a  single-season, 
multiple- state  occupancy  model  (Nichols  et  al.  2007).  We  fit  models  to  the  data  for  three  species: 
Fowler’s  toad  (A naxyrus  fo wleri).  Cope’s  gray  tree  frog  (Hyla  versicolor ),  and  spring  peeper 
( Pseudacris  crucifer).  Before  fitting  the  single-state  model,  we  modified  the  occupancy  data  for 
each  species  by  replacing  all  calling  indices  with  ones,  indicating  simply  that  the  species  was 
detected.  For  the  multiple- state  analyses,  we  treated  the  calling-index  values  as  estimates  of  the 
relative  abundance.  Therefore,  each  wetland  could  be  unoccupied,  occupied  by  a  population 
capable  of  producing  a  low  calling-index  value,  or  occupied  by  a  population  capable  of 
producing  a  high  calling-index  value.  We  pooled  two  values  of  the  calling  index  (either  1  and  2 
or  2  and  3),  which  resulted  in  two  possible  occupancy  states.  We  evaluated  covariates  of 


42 


wetlands  and  the  surrounding  areas  at  three  spatial  extents  that  corresponded  to  breeding, 
migration,  and  dispersal. 

We  evaluated  four  covariates  that  might  be  associated  with  occupancy  and  relative  abundance 
during  the  breeding  phase  of  the  life  cycle:  wetland  area  (m2),  wetland  perimeter,  percentage  of 
canopy  cover  within  a  10-m  buffer  around  the  wetland  perimeter,  and  whether  the  wetland  was 
covered  by  <  50%  or  >  50%  emergent  vegetation  emergent  vegetation. 

We  considered  a  1-km  buffer  around  the  wetland  perimeter  as  the  extent  of  migration.  Within 
that  buffer,  we  estimated  percentage  of  upland  forest  cover,  degree  of  forest  aggregation,  and 
stream  density  (a  proxy  for  riparian  and  mesic-forest  cover)  in  100-m  increments.  We  measured 
forest  aggregation  with  the  clumpiness  index  (FORCLUMP)  in  Fragstats  4.0,  which  estimates 
the  degree  of  fragmentation  of  a  land-cover  type  independent  of  the  area  of  that  land-cover  type. 

We  considered  a  1.5-5  km  buffer  around  the  wetland  perimeter  as  the  extent  of  dispersal.  Within 
that  buffer,  we  estimated  percentage  of  wetland  area,  density  of  highways  with  high  traffic 
volume  (km/km2),  effective  mesh  size  (ha),  and  percentage  of  impervious  surface  in  0.5-km 
increments.  Effective  mesh  size  is  an  estimate  of  the  probability  that  two  random  points  in  a 
given  area  can  be  connected  without  encountering  a  barrier  (in  this  case,  a  road),  and  is 
interpreted  as  the  expected  size  of  a  patch  in  an  area  with  fragmented  land  cover. 

We  first  fit  single-covariate  models  in  which  we  included  values  of  the  covariate  at  all  spatial 
extents  or  increments.  We  considered  an  association  between  the  covariate  and  the  response 
variable  to  be  supported  by  the  data  if  it  was  included  in  a  model  with  an  AICc  value  that  was 
smaller  than  the  intercept-only  model  and  if  the  95%  confidence  interval  around  the  estimates  of 
the  regression  coefficient  for  the  covariate  did  not  include  0.  We  retained  these  covariates  and  fit 
multiple-covariate  models.  We  used  values  of  the  covariate  at  the  distance  and  with  the 
functional  form  that  had  the  lowest  AICc.  We  fit  the  multiple-covariate  models  with  every 
possible  combination  of  two  or  more  covariates.  However,  we  computed  Pearson’s  correlation 
coefficients  for  all  pairs  of  covariates  and  did  not  include  covariates  in  the  same  model  if  the 
absolute  value  of  the  correlation  coefficient  was  >  0.60. 


Responses  of  songbird  occupancy  to  environmental  change 

We  assessed  the  level  of  forest  fragmentation  in  coastal  and  inland  Virginia  on  the  basis  of  30-m 
data  from  the  Landscape  Fire  and  Resource  Management  Planning  Tools  Project  (LANDFIRE) 
(2012).  We  used  the  amount  of  forest  edge  as  a  measure  of  fragmentation.  We  assumed  that 
areas  classified  by  LANDFIRE  as  forest,  grassland,  row  crop,  and  rural  were  relevant  to  deer, 
and  measured  the  mean  length  per  1  km2  plot  of  edges  between  forest  and  grassland,  row  crops, 
or  rural  areas.  We  also  calculated  road  density  (km  per  km2)  in  each  region;  road  density  often  is 
correlated  with  deer  density  (Lovely  et  al.  2013).  As  an  independent  index  of  deer  abundance  for 
the  two  regions,  we  divided  the  number  of  harvested  deer  (2004-2013,  data  provided  by  the 
Virginia  Department  of  Game  and  Inland  Fisheries)  by  forest  area  in  each  county  that  intersected 
our  sites. 
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We  used  two  tests  to  assess  differences  in  deer  use  between  coastal  and  inland  Virginia.  We  used 
a  Wilcoxon  rank  sum  test  to  compare  the  number  of  deer  pellets  in  coastal  and  inland  Virginia, 
and  a  paired  Most  to  compare  the  number  of  deer  harvested.  We  determined  whether  edge 
density  differed  on  the  basis  of  whether  95%  confidence  intervals  overlapped.  We  estimated  both 
bird  and  deer  pellet  densities  (individuals  and  pellets  ha"1,  respectively)  with  Multico variate 
Distance  Sampling  (MCDS)  in  Distance  6.0  software  (Thomas  et  al.  2010).  We  determined 
whether  detection  probabilities  differed  between  species  and  guilds  on  the  basis  of  whether  95% 
confidence  intervals  overlapped.  We  grouped  birds  into  three  guilds:  species  that  nest  in 
understory  shrubs  and  glean  understory  foliage,  species  that  nest  in  the  canopy  and  feed  on  open 
ground,  and  species  that  nest  in  the  canopy  or  cavities  and  either  feed  in  the  canopy  or  sally 
aerially  (Ehrlich  et  al.  1988).  We  restricted  species-level  analyses  to  passerines  that  nest  in  forest 
and  for  which  we  had  >  60  detections  (Thomas  et  al.  2010).  The  latter  restricted  our  species-level 
analysis  to  coastal  Virginia.  We  used  Spearman’s  rank  correlation  to  test  whether  deer  and  bird 
densities  were  correlated. 


Partitioning  drivers  of  beta  diversity 

Calculation  of  turnover  and  nestedness.  We  calculated  beta  diversity  from  detection  data  at 
two  spatial  resolutions:  points  or  transects  and  canyons.  We  calculated  temporal  and  spatial  beta 
diversity  separately.  We  defined  temporal  beta  diversity  as  the  variation  in  species’  identities 
among  years  at  one  point,  transect,  or  canyon.  We  defined  spatial  beta  diversity  as  the  variation 
in  species’  identities  among  all  points  or  transects  within  a  canyon,  or  among  all  canyons  within 
a  mountain  range,  in  a  given  year. 

We  partitioned  beta  diversity  into  two  components — turnover  and  nestedness — with  the  additive 
partitioning  method  outlined  in  Baselga  (2010).  We  based  estimates  of  turnover  and  nestedness 
on  Sprcnsen  dissimilarity,  which  reflects  the  proportion  of  species  shared  between  two 
assemblages.  Additive  partitioning  uses  the  Simpson  dissimilarity  index,  an  estimate  of  beta 
diversity  that  is  insensitive  to  variation  in  species  richness.  Thus,  the  Simpson  dissimilarity  index 
reflects  turnover.  Total  beta  diversity,  defined  as  the  Sprensen  dissimilarity,  comprises  variation 
due  to  turnover  and  variation  due  to  nestedness.  Therefore,  the  difference  between  the  Sprcnsen 
and  Simpson  dissimilarity  indices  is  the  contribution  of  nestedness  to  total  beta  diversity 
(Baselga  2010).  Total  beta  diversity  ranges  from  zero  to  one,  and  turnover  and  nestedness 
together  sum  to  total  beta  diversity.  Therefore,  turnover  and  nestedness  also  range  from  zero  to 
one. 

Because  beta  diversity  estimates  calculated  on  the  basis  of  different  numbers  of  surveys  are  not 
directly  comparable,  we  used  subsampling  to  compare  beta  diversity  estimates  from  different 
locations  or  taxonomic  groups.  For  each  temporal  or  spatial  unit  for  which  we  wished  to 
calculate  beta  diversity  (typically  three  to  ten  surveys),  we  randomly  sampled  100  pairs  of 
surveys  and  used  the  average  pairwise  beta  diversity  over  all  pairs  as  our  estimate  of  beta 
diversity  for  that  set  of  surveys.  The  result  of  this  process  was  estimates  of  temporal  turnover  and 
nestedness  at  each  point  or  transect  and  canyon,  and  estimates  of  spatial  turnover  and  nestedness 
among  all  points  within  each  canyon  and  among  all  canyons  within  each  mountain  range. 
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Null  models  of  turnover  and  nestedness.  We  used  a  null  model  to  estimate  beta  diversity  in  the 
absence  of  any  ecological  structuring  (e.g.,  environmental  filtering  or  limiting  similarity).  For 
each  set  of  surveys  for  which  we  wished  to  calculate  beta  diversity,  we  randomly  sampled  999 
pairs  of  surveys  and,  for  each  pair,  shuffled  the  identities  of  observed  species  while  holding 
species  richness  constant.  We  calculated  beta  diversity  for  each  shuffled  pair  and  used  these  999 
values  to  estimate  an  expected  distribution  of  beta  diversity  values  in  the  absence  of  ecological 
structuring  of  species’  identities. 

We  defined  all  observed  beta  diversity  estimates  less  than  2.5%  or  greater  than  97.5%  of  the 
null-model  estimates  as  evidence  of  ecological  structuring  in  beta  diversity.  Turnover  values  less 
than  2.5%  of  the  null-model  estimates  indicate  more-consistent  species  composition  than  would 
be  expected  in  the  absence  of  ecological  processes,  whereas  values  greater  than  97.5%  of  the 
null-model  estimates  indicate  more-variable  species  composition  than  would  be  expected  in  the 
absence  of  ecological  processes.  Nestedness  estimates  less  than  2.5%  of  the  null-model  estimates 
indicate  less  overlap  in  composition  among  assemblages  than  would  be  expected  in  the  absence 
of  ecological  processes,  whereas  estimates  greater  than  97.5%  of  the  null-model  estimates 
indicate  more  overlap  in  composition  among  assemblages  than  would  be  expected  in  the  absence 
of  ecological  processes. 

Statistical  comparison  of  spatial  resolutions  and  zoogeographic  regions.  We  used  a  Bayesian 
linear  model  to  estimate  differences  in  turnover  and  nestedness  between  the  point  or  transect 
resolution  and  the  canyon  resolution,  and  to  compare  turnover  and  nestedness  between  the 
central  and  western  Great  Basin.  The  model  was 

Vi  =  a  +  |)  vre.st  I(west )  +  P  canyon  I (cany on)  +  Q,  (1) 

where  yn  is  the  observed  turnover  or  nestedness  for  survey  i,  a  is  the  mean  turnover  or  nestedness 
in  the  central  Great  Basin  at  the  point  or  transect  resolution,  P west  is  the  deviation  from  a  in  the 
western  Great  Basin,  I(west )  is  an  indicator  variable  equal  to  one  if  observation  i  is  in  the  western 
Great  Basin  and  zero  otherwise,  ficanyon  is  the  deviation  from  a  at  the  canyon  resolution, 

I(canyon)  is  an  indicator  variable  equal  to  one  if  observation  i  is  estimated  at  the  canyon 
resolution,  and  □  is  the  residual  variation  in  turnover  or  nestedness  for  observation  i. 

We  assigned  unifonn  priors  on  [0,  1]  to  a,  uniform  priors  on  [0  -  a,  1  -  a]  to  ft  west,  and  uniform 
priors  on  [0  -  a  -  fiwest,  1  -  a  -  pM  <est\  to  P canyon.  We  assigned  a  Gaussian  prior  to  Q,  with  zero 
mean  and  standard  deviation  uniformly  distributed  on  [0,  1].  These  priors  reflected  that  turnover 
and  nestedness  estimates  were  between  zero  and  one  in  all  cases.  We  fitted  the  model  defined  in 
eq.  (1)  in  WinBUGS  1.4  with  the  R2WinBUGS  package  in  R.  We  summarized  all  outputs  in  R. 

Statistical  analysis  of  the  relation  between  beta  diversity  and  environmental  variables.  We 

used  a  random  forest  model  to  determine  whether  and  how  turnover  and  nestedness  of  birds  and 
butterflies  in  the  central  and  western  Great  Basin  changed  along  environmental  gradients. 
Random  forests  use  recursive  binary  splits  to  generate  a  large  number  of  tree-like  structures  (a 
forest),  each  of  which  partitions  variation  in  a  response  variable  (e.g.,  turnover  or  nestedness) 
into  groups  on  the  basis  of  a  set  of  covariates  (e.g.,  environmental  variables).  Random  forests 
incorporate  high-level  interactions  and  nonlinear  effects,  so  are  well  suited  to  identifying 
complex  ecological  relations  when  there  are  a  moderate  to  large  number  of  covariates.  The  fitted 
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random  forest  model  is  based  on  the  consensus  output  over  many  trees;  we  used  500  trees  in  all 
analyses. 

We  assessed  model  fit  as  the  r2  value,  based  on  Pearson’s  r,  between  observed  and  fitted  beta 
diversity  estimates.  We  used  ten-fold  cross  validation  to  determine  the  robustness  of  fitted 
relations.  In  cross-validation,  one  partitions  the  data  into  multiple  similar-sized  subsets  of  all 
observations  (folds),  and  fits  the  model  with  each  fold  removed  in  turn  (holdout  data).  Fitted 
models  are  used  to  project  the  response  variable  for  the  holdout  data,  and  cross-validated  r2 
values  are  based  on  the  correlation  between  predicted  and  observed  values.  Cross-validation 
identifies  overfitted  models  because  such  models  generate  poor  predictions  for  holdout  data. 

We  plotted  fitted  models  on  the  basis  of  partial  dependence  plots,  which  estimate  the  marginal 
effect  of  a  given  covariate  on  the  response.  Marginal  effects  are  estimated  by  holding  the  value 
of  the  covariate  constant  while  fitting  random  forest  models  with  a  range  of  values  for  other 
covariates.  The  marginal  effect  is  estimated  as  the  average  association  between  the  covariate  and 
the  response  variable  averaged  over  all  model  runs.  Partial  dependence  plots  do  not  display  the 
multilevel  interactions  supported  by  random  forest  models,  but  approximate  the  shape  of  fitted 
relations.  We  fitted  random  forest  models  with  the  randomForest  package  in  R. 
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Results  and  Discussion 


Hierarchical  estimates  of  species  richness 

The  results  of  this  analysis  provide  information  on  spatial  scales  at  which  it  may  be  most  relevant 
to  assess  and  monitor  native  birds  and  butterflies  and  project  responses  of  these  taxonomic 
groups  to  land  use  and  management,  especially  if  the  response  variable  is  species  richness  or 
occupancy,  throughout  the  Intermountain  West.  Birds  and  butterflies  are  two  of  the  most 
commonly  monitored  faunal  groups,  and  well  over  half  of  the  Intermountain  West  is  managed  by 
federal  agencies,  including  the  Department  of  Defense.  The  results  also  support  the  design  of 
effective,  cost-efficient,  and  spatially  transferable  methods  for  assessment  and  monitoring  of 
native  species;  and  evaluation  of  theories  regarding  relations  between  species  richness  and 
environmental  variables. 

The  central  Great  Basin  bird  data  encompassed  15  years  of  surveys,  with  4-21  points  per  canyon 
and  27  canyons.  The  number  of  years  a  given  point  was  surveyed  ranged  from  6  to  15.  The 
western  Great  Basin  bird  data  encompassed  three  years  of  surveys,  with  8-20  points  per  canyon 
and  13  canyons.  We  treated  canyon  as  an  indicator  variable  and  identified  23  covariates,  nine  at 
the  canyon  level  and  13  at  the  point  level,  that  reasonably  could  be  expected  to  be  associated 
with  species  richness  of  birds.  We  also  included  number  of  survey  years  as  a  covariate. 

At  the  point  level  in  the  central  Great  Basin,  the  relation  between  number  of  survey  years  and 
cumulative  species  richness  was  0.15.  The  adjusted  r-squared  values  of  models  of  point-level 
species  richness  as  a  function  of  point-level  covariates,  canyon  level-covariates,  or  both  were 
0.24,  0.20,  and  0.38,  respectively.  Significant  point-level  covariates  were  elevation  (linear), 
number  of  years  sampled,  canopy  cover,  shrub  cover,  and  mean  minimum  temperature. 
Significant  canyon-level  covariates  were  elevational  range,  terrain  ruggedness  within  a  500-m 
buffer,  canyon  area,  and  riparian  fragmentation.  In  the  model  that  included  covariates  at  both 
spatial  resolutions,  significant  covariates  were  elevation  (linear),  terrain  ruggedness,  mean 
minimum  temperature,  canopy  cover,  years  sampled,  elevational  range  of  the  canyon,  area  of  the 
canyon,  and  riparian  fragmentation.  About  21%  of  the  total  variation  in  cumulative  species 
richness  occurred  among  canyons,  and  about  80%  within  canyons.  The  pseudo-r2  of  the 
hierarchical  model  was  0.41.  The  adjusted  r-squared  of  a  model  of  canyon-level  species  richness 
as  a  function  of  canyon-level  covariates  was  0.86.  Significant  covariates  in  the  latter  model  were 
riparian  fragmentation,  canyon  area,  terrain  ruggedness  within  a  500-m  buffer,  and  sampled  area. 
These  results  suggest  that  birds  may  perceive  and  respond  to  their  environment  at  the  level  of 
canyons  rather  than  at  the  level  of  the  relatively  small  points  typically  used  to  sample  extensive 
areas. 

In  the  western  Great  Basin,  the  adjusted  r-squared  values  of  models  of  point-level  species 
richness  as  a  function  of  point-level  covariates,  canyon  level-covariates,  or  both  were  0.24,  0.14, 
and  0.26,  respectively.  Significant  point-level  covariates  were  elevation  (linear  and  quadratic), 
terrain  ruggedness,  standard  deviation  of  precipitation  and  of  minimum  temperature,  and  canopy 
cover.  Significant  canyon-level  covariates  were  elevational  range,  terrain  ruggedness  within  a 
500-m  buffer,  area  of  canyon  bottom  within  a  500-meter  buffer,  and  number  of  snags.  In  the 
model  that  included  covariates  at  both  spatial  resolutions,  significant  covariates  were  elevation 
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(linear  and  quadratic),  standard  deviation  of  maximum  temperature,  shrub  cover,  elevational 
range  of  the  canyon,  canyon-level  terrain  ruggedness  within  a  500-m  buffer,  area  of  canyon 
bottom  within  a  500-meter  buffer,  and  number  of  snags  in  the  canyon.  About  17%  of  the  total 
variation  in  cumulative  species  richness  occurred  among  canyons,  and  about  83%  within 
canyons.  The  pseudo-r2  of  the  hierarchical  model  was  0.14.  The  adjusted  r-squared  of  a  model  of 
canyon-level  species  richness  as  a  function  of  canyon-level  covariates  was  0.54.  Significant 
covariates  in  the  latter  model  were  elevational  range,  riparian  cover,  and  riparian  fragmentation. 
Again,  canyons  appeared  to  be  a  more  biologically  meaningful  unit  for  breeding  birds  than 
points  within  canyons. 

To  make  an  analysis  of  the  central  Great  Basin  butterfly  data  viable,  we  reduced  the  data  to  two 
years  during  which  the  same  set  of  transects  were  surveyed.  These  reduced  data  included  3-12 
transects  per  canyon  and  15  canyons.  The  western  Great  Basin  butterfly  data  encompassed  4 
years  of  surveys,  with  8-20  transects  per  canyon  and  8  canyons.  We  treated  canyon  as  an 
indicator  variable  and  identified  17  covariates,  six  at  the  canyon  level  and  1 1  at  the  transect  level, 
that  reasonably  could  be  expected  to  be  associated  with  species  richness  of  butterflies. 

In  the  central  Great  Basin,  the  adjusted  r-squared  values  of  models  of  transect-level  species 
richness  as  a  function  of  transect-level  covariates,  canyon  level-covariates,  or  both  were  0.59, 
0.40,  and  0.73,  respectively.  Significant  transect-level  covariates  were  transect  area,  the  standard 
deviation  of  mean  maximum  temperature,  nectar  abundance,  and  mud  abundance.  Significant 
canyon-level  covariates  were  elevational  range,  terrain  ruggedness  within  a  100-m  buffer,  and 
riparian  fragmentation.  In  the  model  that  included  covariates  at  both  spatial  resolutions, 
significant  covariates  were  transect  area,  the  standard  deviation  of  mean  maximum  temperature, 
nectar  abundance,  mud  abundance,  the  area  of  the  canyon  bottom  within  a  500-m  buffer,  and 
terrain  ruggedness  within  a  100-m  buffer.  About  48%  of  the  total  variation  in  cumulative  species 
richness  occurred  among  canyons,  and  about  52%  within  canyons.  Species  richness  was  strongly 
associated  with  transect  area,  nectar  abundance,  mud  abundance,  and  terrain  ruggedness  within  a 
100-m  buffer.  The  pseudo-r2  of  this  hierarchical  model  was  0.60.  We  included  three  canyon- 
level  covariates  in  a  multiple-regression  model  of  canyon-level  species  richness;  the  adjusted  r2 
of  this  model  was  0.57,  and  species  richness  was  most  closely  associated  with  presence  or 
absence  of  riparian  areas. 

In  the  western  Great  Basin,  the  adjusted  r2  values  of  models  of  transect-level  species  richness  as 
a  function  of  transect-level  covariates,  canyon  level-covariates,  or  both  were  0.70,  0.23,  and 
0.71,  respectively.  Significant  transect-level  covariates  were  elevation  (both  linear  and 
quadratic),  mean  precipitation,  mean  minimum  temperature,  nectar  abundance,  and  mud 
abundance.  One  canyon-level  covariate,  elevational  range,  was  statistically  significant.  In  the 
model  that  included  covariates  at  both  spatial  resolutions,  significant  covariates  were  elevation 
(both  linear  and  quadratic),  mean  precipitation,  nectar  abundance,  mud  abundance,  and 
elevational  range  of  the  canyon.  About  22%  of  the  total  variation  in  cumulative  species  richness 
occurred  among  canyons,  and  about  82%  within  canyons.  The  pseudo-r2  of  this  hierarchical 
model  was  0.64.  The  number  of  canyons  was  too  small  for  us  to  model  canyon-level  species 
richness.  However,  our  results  were  consistent  with  the  hypothesis  that  butterflies  perceive  their 
environment  at  smaller  spatial  extents  than  birds. 
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Periodicity  of  sampling:  duration  of  point  counts  of  birds 

The  results  of  this  analysis  provided  knowledge  that  can  facilitate  effective,  cost-efficient,  and 
geographically  transferable  methods  for  sampling  native  birds.  Results  of  analyses  applied  to 
data  from  the  Chesapeake  Bay  Lowlands,  central  Great  Basin,  and  western  Great  Basin  were 
similar  despite  differences  in  ecology  and  land  use  between  the  ecosystems.  The  results  also 
inform  assessment  of  trade-offs  between  increasing  the  duration  of  sampling  and  increasing  the 
number  of  locations  sampled. 

Detection  estimates.  For  82  ±  4%  (mean  ±  SE)  of  the  species  in  a  given  ecosystem  (Chesapeake 
Bay  Lowlands,  central  Great  Basin,  or  western  Great  Basin),  models  of  detection  probability  that 
were  based  on  both  the  5-min  and  8-min  counts  converged  (Table  1).  Detection  probabilities 
based  on  both  the  5-min  and  8-min  counts  were  >  0.3  for  40  ±  4%  of  the  species  in  a  given 
ecosystem  (range  24-54%)  and  <  0.3  for  35  ±  1%  of  the  species  (range  30-40%).  The 
percentage  of  species  with  detection  probabilities  >  0.3  was  higher  in  Chesapeake  Bay  Lowlands 
(54%  in  2012,  47%  in  2013)  than  in  the  central  Great  Basin  (24%  in  2012,  36%  in  2013)  or 
western  Great  Basin  (35%  in  2012,  41%  in  2013).  Increasing  the  count  duration  from  5  min  to  8 
min  increased  the  detection  probability  to  >  0.3  for  6  ±  1%  of  the  species  in  a  given  ecosystem 
(range  3-7%),  but  decreased  the  detection  probability  to  <  0.3  for  2.0  ±1%  (range  =  0-4%)  of 
the  species. 

Occupancy  estimates.  We  estimated  occupancy  for  35  species  in  the  Chesapeake  Bay 
Lowlands,  27  species  in  the  central  Great  Basin,  and  37  species  in  the  western  Great  Basin. 
Occupancy  estimates  based  on  5-min  and  8-min  counts  were  not  statistically  different  for  any  of 
these  species  with  detection  probabilities  >  0.3.  This  suggests  that  a  modest  extension  of  the 
recommended  5-min  duration  of  standardized  counts  (Ralph  et  al.  1993,  1995,  Matsuoka  et  al. 
2014)  is  unlikely  to  affect  inferences  based  on  occupancy  models. 

However,  we  found  that  the  precision  (CV)  of  occupancy  estimates  for  97%  of  the  species  we 
examined  increased  by  12  ±  2%  (range  0-38%)  when  duration  increased  from  5  min  to  8  min. 
This  suggests  a  trade-off  between  count  duration  and  precision  and  that  a  modest  3-min  increase 
in  count  duration  can  provide  more-precise  occupancy  estimates. 

We  could  not  estimate  occupancy  for  46-76%  of  the  species  detected  in  each  ecosystem  and 
year,  and  our  work  highlighted  the  difficulty  of  modeling  occupancy  of  species  with  large  home 
ranges  or  uncommon  species.  Species  with  low  naive  occupancy  (percentage  of  points  occupied, 
not  accounting  for  imperfect  detection)  for  which  models  did  not  converge  included  those  with 
home  ranges  larger  than  the  area  of  the  points  (~3  ha)  (e.g.,  woodpeckers  or  corvids),  those  that 
rarely  occur  in  the  land-cover  types  we  sampled  (e.g.,  synanthropic  and  wetland  species),  and 
those  that  were  rare  in  our  study  areas.  As  a  general  guideline,  MacKenzie  and  Royle  (2005) 
recommended  increasing  the  number  of  counts  rather  than  the  number  of  points  to  maximize  the 
precision  of  estimates  for  common  species,  and  increasing  the  number  of  points  rather  than  the 
number  of  counts  to  maximize  precision  for  rare  species. 
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Table  1.  Total  number  of  species  detected  in  2012  and  2013  in  the  Chesapeake  Lowlands,  central  Great  Basin,  and  western  Great 
Basin;  the  number  (percentage)  of  species  for  which  models  did  or  did  not  converge;  and  the  number  (percentage)  of  species  for  which 
detection  probabilities  (p)  based  on  5-min  and  8-min  counts  were  either  >  0.3  or  <  0.3.  Percentages  are  based  on  the  total  number  of 
species  detected  within  a  given  year  and  ecosystem. 


Chesapeake  Bay 
Lowlands 

central  Great 
Basin 

western  Great 
Basin 

2012 

2013 

2012 

2013 

2012 

2013 

Number  of  species  detected 

57 

59 

70 

73 

75 

76 

Models  based  on  5-min  and  8-min  counts  converged 

53 (93) 

54  (92) 

51  (73) 

57  (78) 

55 (73) 

61  (80) 

Models  converged  when  based  on  8-min  but  not  5-min 
counts 

0(0) 

3(5) 

1(1) 

5(7) 

1(1) 

0(0) 

Models  converged  when  based  on  5-min  but  not  8-min 
counts 

4(7) 

2(3) 

1(1) 

0(0) 

2(3) 

0(0) 

Neither  model  converged 

0(0) 

0(0) 

17  (24) 

11(15) 

17  (23) 

15  (20) 

p  >  0.3  when  based  on  5-min  and  8-min  counts 

31  (54) 

28  (47)a 

17  (24) 

26  (36) 

26  (35) 

31  (41) 

p  >  0.3  when  based  on  8-min  but  not  5-min  counts 

4(7) 

4(7) 

4(6) 

4(5) 

2(3) 

4(5) 

p  >  0.3  when  based  on  5-min  but  not  8-min  counts 

1(2) 

1(2) 

2(3) 

3(4) 

1(1) 

0(0) 

p  <  0.3  when  based  on  5-min  and  8-min  counts 

17  (30) 

21  (36) 

28  (40) 

24  (33) 

26  (35) 

26  (34) 

a  Blue-gray  Gnatcatchers  excluded  because  confidence  intervals  for  occupancy  ranged  from  zero  to  one. 
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Across  ecosystems  and  years,  point  estimates  of  occupancy  for  19  ±  5%  of  the  species  in  a  given 
ecosystem  were  0.001-0.32  higher  when  based  on  5-min  than  8-min  counts.  This 
counterintuitive  result  reflects  the  manner  in  which  occupancy  is  estimated.  Occupancy  estimates 
based  on  5-min  counts  can  be  higher  than  those  based  on  8-min  counts  when  naive  occupancy 
from  5-min  and  8-min  counts  are  similar,  but  the  species  is  detected  on  fewer  of  the  5-min  counts 
than  the  8-min  counts.  Although  the  difference  in  occupancy  estimates  was  slight,  we  suspect 
that  it  might  be  sufficient  to  affect  estimates  of  species  richness  based  on  occupancy  models 
(Iknayan  et  al.  2014),  especially  because  sample  sizes  were  too  small  to  estimate  occupancy  for  a 
substantial  proportion  of  the  avifauna  in  our  ecosystems. 

Hayes  and  Monfils  (2015)  suggested  that  point-count  data  are  not  well  suited  for  occupancy 
modeling  because  the  mobility  of  birds  violates  the  closure  assumption — the  occupancy  status  of 
points  does  not  change  among  counts.  Because  the  area  of  most  points  is  smaller  than  the 
smallest  home  range  of  the  species  sampled,  individuals  regularly  move  in  and  out  of  the  points 
(Rota  et  al.  2009).  Simulations  indicated  that  even  slight  movements  resulted  in  detection 
probabilities  that  were  lower  than  the  true  value  and  occupancy  estimates  that  were  higher  than 
the  true  value  (Hayes  and  Monfils  2015).  Our  results  were  consistent  with  the  results  of  these 
simulations.  Ultimately,  whether  point-count  data  can  be  used  to  estimate  occupancy  hinges  on 
the  definition  of  use  (i.e.,  open  system)  versus  occupancy  (i.e.,  closed  system)  (Latif  et  al.  2016). 


Periodicity  of  sampling:  trade-offs  between  single-day  and  multiple-days  sampling 

The  results  of  this  analysis  provide  knowledge  that  can  facilitate  effective,  cost-efficient,  and 
geographically  transferable  methods  for  sampling  native  anurans,  birds,  and  butterflies. 
Inferences  about  trade-offs  between  single-day  and  multiple-days  sampling  generally  were 
consistent  among  taxonomic  groups  and  geographic  regions  despite  substantial  differences  in 
biology  among  the  taxonomic  groups  and  ecology  and  land  use  among  the  regions.  Additionally, 
these  analyses  provide  evidence  for  evaluating  whether  the  assumptions  of  occupancy  theory  are 
met  in  diverse  taxonomic  groups  and  ecosystems. 

Birds.  After  removing  species  that  were  not  detected  during  one  of  the  two  years,  were  not 
detected  in  a  sufficient  number  of  points,  or  for  which  models  of  detection  probability  did  not 
converge,  we  were  able  to  analyze  data  for  23  species  in  2012  and  22  species  in  2013. 

Occupancy  estimates  derived  from  the  multiple-days  design  were  higher  than  those  derived  from 
the  single-day  design  for  all  species  in  2012  (mean  0.17,  range  0.04-0.33)  and  2013  (mean  0.16, 
range  0.01-0.38).  However,  95%  confidence  intervals  from  the  multiple-days  and  single-day 
models  overlapped  for  20  species  in  2012  and  all  species  in  2013.  Estimates  of  p*  were  higher 
when  based  on  data  from  the  single-day  design  than  from  the  multiple-days  design  for  all  species 
in  2012  (mean  0.14,  range  0.01-0.39)  and  21  species  in  2013  (mean  0.1 1,  range  -  0.09-0.33). 
However,  95%  confidence  intervals  around  estimates  of  p*  based  on  the  two  sampling  designs 
overlapped  for  9  species  in  2012  and  19  species  in  2013. 

There  are  at  least  two  reasons  why  occupancy  estimates  based  on  the  multiple-days  design 
consistently  were  higher  than  those  based  on  the  single-day  design.  First,  multiple-days  estimates 
of  detection  probability  were  lower  than  the  single-day  estimates.  As  detection  probability 
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decreases,  detection-weighted  occupancy  estimates  increase  relative  to  naive  occupancy 
estimates.  Second,  the  proportion  of  sites  on  which  species  were  detected  during  at  least  one 
survey  was  higher  when  the  multiple-days  design  was  implemented  than  when  the  single-day 
design  was  implemented.  These  two  differences  in  detection  probability  between  the  two 
sampling  designs  likely  would  be  consistent  in  other  regions  and  time  periods. 

Our  results  suggest  that  the  closure  assumption  was  violated  in  the  multiple-days  design.  Our 
findings  are  consistent  with  those  of  Rota  et  al.  (2009),  who  found  that  the  closure  assumption 
was  violated  for  all  1 8  species  of  breeding  songbirds  that  they  surveyed  over  8  days  in  deciduous 
forests  in  the  eastern  United  States.  The  multiple-days  and  single-day  designs  may  estimate 
different  parameters  in  species  with  a  high  proportion  of  mobile  individuals.  Because  the 
multiple-days  design  is  more  likely  than  the  single-day  design  to  sample  both  breeding  and 
dispersing  males,  we  suggest  that  it  estimates  the  probability  of  use  rather  than  of  occupancy. 

The  period  of  sampling  in  the  single-day  design  is  short,  which  suggests  that  this  design  has  a 
relatively  high  probability  of  meeting  the  closure  assumption.  However,  the  single-day  design  is 
less  likely  than  the  multiple-days  design  to  sample  dispersing  individuals  or  individuals  that  are 
not  associated  with  a  territory. 

Butterflies.  Analyses  are  ongoing. 

Anurans.  After  removing  species  that  did  not  meet  our  criteria  for  consistent  detection,  we  were 
able  to  analyze  two  years  of  data  for  5  species  that  breed  early  in  the  season,  7  species  that  breed 
during  the  middle  of  the  season,  and  4  species  that  breed  late  in  the  season.  Some  individual 
species  breed  over  a  relatively  long  period  of  time,  and  therefore  were  included  in  more  than  one 
analysis.  Occupancy  estimates  derived  from  the  multiple-days  design  were  higher  than  those 
derived  from  the  single-day  design  for  all  16  species  models  in  2012  (mean  0.14,  range  0.02- 
0.47)  and  13  of  14  species  in  2014  (mean  0.12,  range  =  -0.007-0.35).  However,  95%  confidence 
intervals  around  estimates  of  occupancy  overlapped  in  all  cases. 

Estimates  of  p*  were  higher  when  based  on  data  from  the  single-day  design  than  from  the 
multiple-days  design  in  all  cases  (2012:  mean  0.09,  range  =  0.0002-0.28;  2014:  mean  0.14,  range 
0.002-0.66).  However,  95%  confidence  intervals  around  estimates  of  p*  based  on  the  two 
sampling  designs  overlapped  for  7  species  in  2012  and  8  species  models  in  2014. 

Simulations.  The  extent  to  which  availability  affected  estimates  of  occupancy  and  p*  differed 
between  survey  designs  (Figure  9).  Availability  led  to  greater  bias  in  occupancy  estimates  in  the 
single-day  than  in  the  multiple-days  design.  Occupancy  estimates  based  on  the  multiple-days 
were  unbiased  when  availability  was  >  0.4,  but  occupancy  estimates  based  on  the  single-day 
design  were  biased  whenever  availability  was  <  1.0.  By  contrast,  availability  led  to  greater  bias 
in  estimates  of  p *  in  the  multiple-days  design  than  in  the  single-day  design.  Estimates  of  p* 
based  on  the  single-day  design  were  unbiased  when  availability  was  >  0.2,  but  estimates  of  p* 
based  on  the  multiple-days  design  were  biased  whenever  availability  was  <  1.0. 
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Figure  9.  Simulation-based  responses  of  estimates  of  occupancy  and  probability  of  detecting  a 
given  species  at  least  once  (p*)  to  changes  in  availability  of  the  species. 

Although  our  work  suggests  that  the  closure  assumption  is  met  in  the  single-day  design,  this 
design  potentially  violates  two  other  assumptions  of  occupancy  models.  First,  it  may  violate  the 
assumption  of  independent  observations  because  an  observer  tends  to  remember  species  detected 
during  previous  surveys,  particularly  if  the  same  observer  conducts  all  surveys  (MacKenzie  et  al 
2006).  Nevertheless,  for  nearly  all  species  in  our  analyses,  models  that  included  an  effect  of 
previous  detections  on  detection  probability  were  not  strongly  supported.  Second,  the  single-day 
design  may  introduce  heterogeneity  in  detection  probabilities  among  sites  (MacKenzie  et  al. 
2006)  because  different  groups  of  sites  are  surveyed  on  different  days,  and  environmental 
conditions  may  vary  among  days. 

The  precision  of  estimates  of  occupancy  and  detection  probability  that  were  based  on  the  single¬ 
day  design  was  higher  than  those  based  on  the  multiple-days  design.  Higher  precision  reflected  a 
lower  proportions  of  locations  at  which  a  species  was  observed  during  a  subset  of  surveys. 
Increasing  the  number  of  surveys  at  a  given  site  typically  increases  precision  (MacKenzie  and 
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Royle  2005),  but  there  is  a  trade-off  between  the  number  of  surveys  conducted  at  a  given  site  and 
the  number  of  sites.  In  many  ecosystems,  the  single-day  design  allows  greater  flexibility  in 
allocation  of  the  number  of  sites  and  surveys.  All  else  being  equal,  precision  of  occupancy  and 
detection  estimates  for  common  species  will  be  maximized  by  conducting  a  greater  number  of 
surveys  at  fewer  sites,  and  precision  of  estimates  for  rare  species  will  be  maximized  by 
conducting  fewer  surveys  at  a  greater  number  of  sites  (MacKenzie  and  Royle  2005). 

Designing  sampling  to  estimate  occupancy  requires  consideration  of  the  sampling  effort 
available  and  trade-offs  in  allocation  of  that  effort  between  numbers  of  sites  and  surveys 
(MacKenzie  and  Royle  2005).  In  most  cases,  less  time  is  spent  traveling  between  sites  on  a  given 
day  when  the  single-day  design  is  implemented  than  when  the  multiple-days  design  is 
implemented  (MacKenzie  et  al.  2006,  Lele  2012).  Because  breeding  birds  cannot  be  sampled  for 
more  than  four  or  five  hours  after  sunrise,  the  single-day  design  generally  allows  one  to  sample  a 
greater  number  of  sites  than  the  multiple-days  design.  We  found  that  conducting  three  surveys  of 
131  sites  within  14  parks  with  a  multiple-days  design  required  approximately  twice  as  much  time 
(214  hr)  as  did  sampling  with  the  single-day  design  (1 10  hr).  These  estimates  did  not  include 
return-travel  time,  which  did  not  limit  sampling  effort.  As  a  result,  we  could  visit  more  sites  or 
conduct  a  greater  number  of  surveys  per  site  by  implementing  the  single-day  design.  If  the 
objective  is  to  define  habitat  quality  for  birds,  then  we  suggest  implementing  the  multiple-days 
design  because  the  estimates  of  occupancy  derived  from  this  design  reflect  habitat  use 
(MacKenzie  et  al.  2006).  If  the  objective  is  to  monitor  site  or  patch  occupancy  of  birds,  then  we 
suggest  implementing  the  single-day  design  because  per-site  sampling  effort  is  lower  and 
precision  of  occupancy  estimates  is  higher  than  when  the  multiple- visit  design  is  applied. 


Estimation  of  species  richness  of  multiple  taxonomic  groups  on  the  basis  of  a  single  group 

The  results  of  this  analysis  provide  knowledge  that  can  facilitate  effective,  cost-efficient,  and 
geographically  transferable  methods  for  assessing  and  monitoring  species  richness  of  native 
faunas  and  projecting  responses  of  these  taxonomic  groups  to  land  use  and  management.  The 
methods  that  we  developed  can  be  applied  to  any  taxonomic  group  in  any  location,  although  the 
extent  to  which  indicator  species  are  likely  to  be  effective  well  may  vary  among  taxonomic 
groups  and  locations.  Additionally,  these  analyses  provide  evidence  for  evaluating  theory  on 
whether  objectively  selected  subsets  of  a  given  taxonomic  group  support  inference  to  species 
richness,  whether  currently  or  in  other  time  periods  or  locations. 

Our  exploration  of  potential  surrogates  for  estimation  of  species  richness  is  novel  in  the  extent  to 
which  we  externally  evaluated  both  spatial  and  temporal  transferability.  Our  models  were  robust 
to  the  spatial  and  temporal  variation  in  sampling  locations  that  is  relatively  common  in  long-term 
data.  Furthermore,  we  found  that  the  predictive  accuracy  of  models  based  on  indicator  species 
could  exceed  that  of  models  based  on  environmental  variables. 

Point-level  species  richness  of  birds  ranged  from  1-22,  and  mean  point-level  species  richness 
was  greater  in  the  western  than  in  the  central  Great  Basin  (12.0  ±  3.5  [mean  ±  SD]  and  8.2  ±  3.6, 
respectively).  Similarly,  mean  annual  species  richness  was  greater  in  the  western  Great  Basin 
(89.3  ±  4.9  from  2012-2014)  than  in  the  central  Great  Basin  (81.0  ±  2.6  from  2012-2014).  Mean 
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species  richness  of  butterflies  at  a  given  transect  in  the  central  Great  Basin  was  about  twice  that 
in  the  western  Great  Basin  (21.6  ±  9.5  and  11.5  ±  7.2,  respectively),  but  mean  annual  species 
richness  of  butterflies  was  greater  in  the  western  (e.g.,  82.3  ±  4.0  from  2013-2014)  than  in  the 
central  Great  Basin  (72.5  ±  3.5). 

Explanatory  ability  of  fitted  models.  Before  confronting  our  models  with  independent  data — 
data  from  locations  or  time  periods  that  were  not  used  to  build  the  models — we  evaluated  the 
extent  to  which  the  models  explained  variation  in  observed  species  richness.  Explained  variation 
does  not  provide  information  on  whether  the  model  is  likely  to  accurately  predict  species 
richness  in  those  new  locations  or  time  periods.  Fitted  values  of  species  richness  that  were  based 
on  both  indicator  species  and  environmental  variables  explained  0.27-0.69  of  the  variation  in 
observed  species  richness  of  birds  and  0.52-0.78  of  the  variation  in  observed  species  richness  of 
butterflies  (Table  2).  Models  that  were  based  on  same-taxon  indicator  species  explained  more 
variation  than  models  that  included  cross-taxonomic  indicator  species.  Apart  from  one  case  in 
which  the  explained  variation  was  essentially  equal,  models  that  were  based  on  both  species  and 
environmental  variables  explained  0.07-0.40  more  variation  in  species  richness  than  species- 
only  models,  and  0-0.23  more  variation  than  models  that  were  based  on  environmental  variables 
only. 

The  accuracy  of  these  predictions  may  be  considerable  statistically  (Mpller  and  Jennions  2002), 
but  whether  the  predictions  are  sufficient  to  inform  decision  making  depends  on  the  particular 
situation  (Runge  et  al.  2011,  Keisler  et  al.  2013).  It  is  possible  that  local  managers  or  citizens 
more  readily  could  be  trained  to  sample  a  relatively  small  number  of  species  than  to  sample  all 
species  in  a  fairly  large  taxonomic  group.  We  strongly  recommend  external  evaluation  of  any 
models  that  might  be  used  to  inform  decision-making. 
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Table  2.  Variation  in  species  richness  of  birds  or  butterflies  correctly  explained  by  fitted  random  forest  models  that  were  based  on 
indicator  species  and  environmental  variables  (i.e.,  the  fit  of  the  random  forest  model  to  the  data  used  for  its  construction).  Central  and 
western  refer  to  subregions  of  the  Great  Basin,  r2  values  are  based  on  Pearson’s  r  and  were  averaged  across  the  ten  iterations  of  each 
model.  Full  models  included  species  and  environmental  variables. 


Taxonomic 
group  of 
response 
variable 

Taxonomic  group 

of  indicator  species  Subregion  and  years 

r2,  full  model 

r2,  species-only 
model 

r2,  environmental 
variables-only  model 

birds 

birds 

western  2012-2014 

0.50 

0.32 

0.43 

birds 

butterflies 

western  2012-2014 

0.39 

0.09 

0.39 

birds 

birds 

central  2008-2010 

0.57 

0.48 

0.46 

birds 

butterflies 

central  2001-2003 

0.27 

0.20 

0.07 

birds 

birds 

central  2008-2014 

0.69 

0.70 

0.56 

birds 

butterflies 

central  2013-2014 

0.55 

0.15 

0.55 

birds 

birds 

western  2012 

0.44 

0.37 

0.33 

birds 

butterflies 

western  2012 

0.30 

0.07 

0.30 

butterflies 

butterflies 

western  2012-2015 

0.78 

0.64 

0.68 

butterflies 

birds 

western  2012-2014 

0.73 

0.38 

0.73 

butterflies 

butterflies 

central  1995-1997 

0.78 

0.67 

0.55 

butterflies 

birds 

central  2001-2003 

0.52 

0.18 

0.42 

butterflies 

butterflies 

central  2001-2014 

0.71 

0.53 

0.61 

butterflies 

birds 

central  1995-2014 

0.54 

0.31 

0.49 

butterflies 

butterflies 

western  2012-2013 

0.73 

0.56 

0.62 

butterflies 

birds 

western  2012 

0.60 

0.41 

0.52 
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Fitted  associations  between  species  richness  and  environmental  variables.  Species  richness 
of  birds  usually  increased  as  canopy  cover  and  the  proportion  of  riparian  cover  increased,  and 
decreased  as  riparian  fragmentation  increased  (Table  3).  The  association  between  species 
richness  of  birds  and  canyon  area  was  unimodal  and  concave  (U-shaped).  The  shape  of 
associations  between  species  richness  of  birds  and  number  of  trees,  terrain  ruggedness,  and 
minimum  temperature  varied  noticeably  between  the  central  and  western  Great  Basin.  Species 
richness  of  butterflies  generally  increased  as  the  maximum  temperature  decreased  and  as  the 
abundance  of  nectar  and  mud  increased.  Species  richness  of  butterflies  also  increased  as  the 
sampled  area  increased,  and  appeared  to  reach  an  asymptote  in  the  western  Great  Basin.  In  the 
central  Great  Basin,  species  richness  of  butterflies  increased  as  mean  elevation  decreased  and  as 
precipitation  increased.  In  the  western  Great  Basin,  by  contrast,  species  richness  was  greatest  at 
intermediate  values  of  elevation  and  precipitation. 

Model  evaluation.  Consistent  with  our  previous  work,  more  variation  in  species  richness  of  a 
given  taxonomic  group  was  predicted  by  same-taxon  indicator  species  than  by  cross-taxonomic 
indicator  species  (Thomson  et  al.  2007).  Much  more  variation  in  species  richness  of  birds  was 
predicted,  whether  in  space  or  in  time,  by  same-taxon  indicator  species  than  by  cross-taxonomic 
indicator  species  (by  an  additional  0.12-0.40  in  full  models,  and  0.32-0.51  in  species-only 
models).  The  most  accurate  model,  which  was  built  with  data  from  the  central  Great  Basin  and 
tested  against  data  from  a  later  time  period  in  the  same  region,  predicted  0.58  of  the  variation  in 
species  richness  of  birds.  Whether  the  full  model  or  species-only  model  predicted  more  variation 
in  species  richness  was  inconsistent,  but  differences  between  the  two  model  types  generally 
increased  as  the  total  variation  explained  decreased.  Two  of  the  environmental  variables-only 
models  predicted  >0.25  of  the  variation  in  species  richness  of  birds. 

With  one  exception,  same-taxon  indicator  species  predicted  more  variation  in  species  richness  of 
butterflies  than  cross-taxonomic  indicator  species  (0.15-0.61  more  in  full  models,  0.08-0.68  in 
species-only  models)  (Table  2).  However,  occurrence  patterns  of  birds  predicted  a  maximum  of 
0.65  of  the  variation  in  species  richness  of  butterflies,  whereas  occurrence  patterns  of  butterflies 
predicted  no  more  than  0.19  of  the  variation  in  species  richness  of  birds.  The  most  accurate 
model,  which  predicted  0.67  of  the  variation  in  species  richness  of  butterflies  in  the  western 
Great  Basin  in  2014,  was  built  with  data  on  environmental  variables  (Table  2).  Four  other 
environmental  variables-only  models  predicted  >0.25  of  the  variation  in  species  richness  of 
butterflies. 

A  possible  explanation  for  the  relatively  poor  predictions  of  species  richness  of  birds  from 
butterfly  indicator  species  is  that  detection  probabilities  and  occupancy  of  adult  butterflies  often 
increase  as  the  abundance  of  nectar  or  mud  increases  (Fleishman  et  al.  2017),  and  in  the  Great 
Basin,  these  resources  generally  are  most  abundant  in  riparian  areas  with  moderate  canopy  cover. 
Relatively  few  butterflies  in  the  Great  Basin  occur  in  comparatively  xeric  shrublands  or 
coniferous  woodlands,  which  provide  habitat  for  numerous  species  of  birds. 

In  external  evaluation,  models  based  on  same-taxon  indicator  species  predicted  more  variation  in 
species  richness  than  those  based  on  environmental  variables  in  seven  of  eight  cases. 
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Table  3.  Variation  in  species  richness  of  birds  or  butterflies  correctly  predicted  by  indicator  species  and  environmental  variables. 
Central  and  western  refer  to  subregions  of  the  Great  Basin  (see  text  and  Fig.  1).  r2  values  are  based  on  Pearson’s  r.  Full  models 
included  species  and  environmental  variables.  Grey  shading  indicates  models  that  predicted  <0.25  of  the  variation  in  species  richness. 


Taxonomic 
group  of 
response 
variable 

Taxonomic 
group  of 
indicator 
species 

Subregion  and  years 
used  to  build  the 
model 

Subregion  and  years 
used  to  test  predictions 
of  the  model 

r2,  full  model 

r2,  species- 
only  model 

r2,  environmental 
variables-only 
model 

birds 

birds 

western  2012-2014 

central  2012-2014 

0.41 

0.47 

0.26 

birds 

butterflies 

western  2012-2014 

central  2013-2014 

0.18 

0.02 

0.17 

birds 

birds 

central  2008-2010 

central  2012-2014 

0.58 

0.52 

0.46 

birds 

butterflies 

central  2001-2003 

central  2013-2014 

0.19 

0.01 

0.20 

birds 

birds 

central  2008-2014 

western  2012-2014 

0.46 

0.50 

0.11 

birds 

butterflies 

central  2013-2014 

western  2012-2014 

0.06 

0.01 

0.07 

birds 

birds 

western  2012 

western  2014 

0.30 

0.34 

0.23 

birds 

butterflies 

western  2012 

western  2014 

0.18 

0.02 

0.17 

butterflies 

butterflies 

western  2012-2015 

central  2013-2014 

0.59 

0.60 

0.30 

butterflies 

birds 

western  2012-2014 

central  2013-2014 

0.11 

0.01 

0.23 

butterflies 

butterflies 

central  1995-1997 

central  2013-2014 

0.70 

0.70 

0.12 

butterflies 

birds 

central  2001-2003 

central  2013-2014 

0.09 

0.02 

0.03 

butterflies 

butterflies 

central  2001-2014 

western  2012-2014 

0.51 

0.57 

0.27 

butterflies 

birds 

central  1995-2014 

western  2012-2014 

0.36 

0.01 

0.39 

butterflies 

butterflies 

western  2012-2013 

western  2015 

0.53 

0.48 

0.46 

butterflies 

birds 

western  2012 

western  2014 

0.65 

0.40 

0.67 
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Environmental  variables-only  models  predicted  0.01-0.38  more  variation  than  species-only 
models  in  the  seven  cases  in  which  the  latter  predicted  <0.02  variation  in  species  richness.  Thus, 
although  there  were  instances  in  which  the  proportion  of  variation  predicted  by  environmental 
variables  was  greater  than  the  proportion  predicted  by  indicator  species,  most  of  these  models 
predicted  relatively  little  variation  in  species  richness  in  an  absolute  sense. 

Identity  and  predictive  capacity  of  indicator  species.  Thirty-four  species  of  birds  were 
included  in  >1  iteration  of  at  least  one  of  the  four  spatially  transferable  models  of  the  species 
richness  of  birds  (i.e.,  two  full  and  two  indicator-species  only  models).  However,  there  was  some 
consistency  in  the  identity  of  the  species  included  in  multiple  models.  Dusky  Flycatcher 
( Empidonax  oberholseri).  Mountain  Chickadee  ( Poecile  gambeli),  Bushtit  ( Psaltriparus 
minimus ),  House  Wren  ( Troglodytes  aedon),  American  Robin  ( Turdus  migratorius ),  Western 
Tanager  (Pi  rang a  ludoviciana),  and  Cassin’s  Finch  ( Haemorhous  cassinii )  were  included  in  >1 
iteration  of  all  of  the  spatially  transferable  models. 

Thirty  species  of  birds  were  included  in  >1  iteration  of  at  least  one  of  the  four  temporally 
transferable  models  of  the  species  richness  of  birds.  The  indicator  species  included  in  the  greatest 
proportion  of  models  of  the  species  richness  of  birds,  all  of  which  had  positive  associations  with 
species  richness,  are  fairly  common  and  widespread  in  both  the  central  and  western  Great  Basin. 
Dusky  Flycatchers  are  associated  with  riparian  areas  and  adjacent  vegetation;  House  Wrens  with 
woodlands,  often  near  edges;  American  Robins  with  woodland  edges;  and  Cassin’s  Finches  in 
open,  coniferous  woodlands  and  in  sagebrush  ( Artemisia  spp.)  steppe  with  scattered  trees,  often 
at  relatively  high  elevations.  Our  models  predict  that  the  locations  in  which  these  four  species  co¬ 
occur  (most  of  which  are  riparian)  have  relatively  high  species  richness.  Sage  Thrasher  and 
Sagebrush  Sparrow  ( Artemisiospiza  nevadensis),  which  had  strong  negative  associations  with 
species  richness  of  birds,  are  associated  with  large  sagebrush  in  areas  with  little  tree  cover.  These 
species  are  most  common  in  the  valleys  separating  mountain  ranges;  we  sampled  relatively  few 
locations  in  which  they  were  likely  to  occur,  and  recorded  few  individuals. 

Twenty-six  species  of  butterflies  were  included  in  >1  iteration  of  at  least  one  of  the  four  spatially 
transferable  models  of  the  species  richness  of  butterflies.  As  with  birds,  there  was  some 
consistency  in  the  identity  of  species  included  in  multiple  models.  Twenty-four  species  of 
butterflies  were  included  in  >1  iteration  of  at  least  one  of  the  four  temporally  transferable  models 
of  the  species  richness  of  butterflies.  Both  of  the  butterfly  species  that  were  included  in  a  high 
proportion  of  models  of  the  species  richness  of  butterflies  are  fairly  widespread  and  conspicuous. 
Icaricia  lupini  occurs  in  diverse  land-cover  types.  Papilio  rutulus  generally  is  associated  with 
riparian  areas,  although  it  is  quite  vagile.  I.  lupini  has  multiple  generations  per  season,  whereas 
P.  rutulus  is  univoltine  and  generally  flies  during  the  middle  of  the  season. 

Thirty- six  species  of  birds  were  included  in  at  least  one  of  the  two  spatially  and  two  temporally 
transferable  models  of  the  species  richness  of  butterflies  that  explained  >0.25  of  the  variation  in 
the  response  variable.  The  bird  species  that  had  the  strongest  negative  associations  with  species 
richness  of  butterflies  typically  occur  in  relatively  dense  woodlands  or  in  shrub-dominated 
uplands,  which  generally  are  relatively  late-successional  and  xeric  communities  with  few  sources 
of  nectar.  However,  not  all  species  with  such  land-cover  associations  had  negative  relations  with 
species  richness  of  butterflies.  For  example,  Sage  Thrasher  ( Oreoscoptes  montanus ),  which  had  a 
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fairly  strong  positive  association  with  species  richness  of  butterflies,  occurs  in  shrub-dominated 
uplands. 

Strengths  of  association  between  species  richness  and  the  occurrence  of  butterfly  indicator 
species  generally  were  appreciably  higher  than  those  between  species  richness  and  the 
occurrence  of  bird  indicator  species. 

Our  work  may  suggest  that  tracking  the  occurrence  of  a  subset  of  an  assemblage  during  periods 
of  environmental  change  will  allow  inference  to  species  richness  of  the  assemblage.  However, 
given  that  adaptive  phenotypic  plasticity  and  evolution  are  species-specific  (Ghalambor  et  al. 
2007,  Kinnison  and  Hariston  2007,  Reed  et  al.  201 1),  current  relations  between  indicator  species 
and  species  richness  likely  are  not  stationary. 


Use  of  occupancy  as  a  surrogate  measure  of  species  richness:  hierarchical  models  of 
occupancy 

Analyses  are  ongoing.  Results  of  the  analyses  will  provide  knowledge  on  the  extent  to  which 
groups  of  birds  with  similar  nesting  habitat  are  likely  to  differ  in  their  response  to  environmental 
change,  and  whether  these  differences  are  consistent  among  ecosystems.  Results  also  will  inform 
selection  of  environmental  variables  for  inclusion  in  assessment  and  monitoring  of  species 
richness. 

Preliminary  results  suggest  that  in  the  Chesapeake  Bay  Lowlands,  occupancy  of  species  that  nest 
on  the  ground  or  in  low  shrubs  is  negatively  associated  with  number  of  snag  saplings  and  mean 
canopy  height,  whereas  occupancy  of  the  same  group  is  positively  associated  with  number  of 
deciduous  saplings  and  proportion  of  riparian  cover.  Persistence  of  species  that  nest  on  the 
ground  or  in  low  shrubs  appears  to  be  positively  associated  with  proportion  of  mesic  deciduous 
or  mixed  tree  cover  and  proportion  of  conifer  cover.  Occupancy  of  species  that  nest  in  tall  trees 
or  shrubs  was  positively  associated  with  structural  heterogeneity,  and  persistence  of  this  group 
was  positively  associated  with  proportion  of  developed  land.  The  latter  may  reflect  planting  of 
large  ornamental  trees  in  many  suburban  areas.  Occupancy  of  species  that  nest  in  tree  cavities 
seems  to  be  negatively  associated  with  number  of  deciduous  saplings  (which  are  too  small  to 
support  construction  of  a  canopy  nest)  and  positively  associated  with  proportion  of  riparian 
cover. 

Preliminary  results  suggest  that  detection  probability  of  all  groups  in  the  western  Great  Basin  is 
negatively  associated  with  terrain  ruggedness.  Complex  terrain  may  interfere  with  both  spatial 
transmission  of  vocalizations  and  visual  detection  of  birds.  Occupancy  of  species  that  nest  on  the 
ground  or  in  low  shrubs  may  be  negatively  associated  with  riparian  fragmentation  and  positively 
associated  with  elevation.  Persistence  of  species  that  nest  on  the  ground  or  in  low  shrubs  seems 
to  be  lower  in  the  Sierra  Nevada  than  the  Wassuk  or  Sweetwater  Range,  and  colonization  seems 
to  be  negatively  associated  with  canopy  prevalence.  Occupancy  of  species  that  nest  in  tall  trees 
or  shrubs  was  positively  associated  with  canopy  prevalence.  Persistence  of  the  latter  group  was 
negatively  associated  with  riparian  fragmentation  and  positively  associated  with  riparian  cover 
and  canopy  prevalence.  Occupancy  of  species  that  nest  in  tree  cavities  was  negatively  associated 
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with  elevation  and  positively  associated  with  terrain  ruggedness,  whereas  persistence  of  the 
group  was  positively  associated  with  canopy  prevalence  and  elevation.  Results  from  the  central 
Great  Basin  are  forthcoming. 


Use  of  occupancy  as  a  surrogate  measure  of  species  richness:  estimation  of  the  occupancy 
of  butterflies  in  diverse  biogeographic  regions 

Results  of  these  analyses  provide  knowledge  on  the  extent  to  which  different  species  of 
butterflies  may  differ  in  their  response  to  environmental  change,  and  whether  these  differences 
are  consistent  among  ecosystems.  Our  results  also  provide  information  to  evaluate  assumptions 
of  occupancy  theory. 

Geographic  variation  in  the  extent  to  which  assumptions  of  occupancy  theory  are  met,  and  the 
manner  in  which  such  variation  may  be  associated  with  regional  differences  in  the  biology  of  the 
focal  species  or  taxonomic  group,  generally  has  not  been  addressed  in  the  literature.  Pursuing 
these  areas  of  inquiry  can  inform  the  design  of  effective  and  cost-efficient  assessment  and 
monitoring.  We  found  that  the  extent  to  which  the  closure  assumption  was  met  varied  among 
ecosystems  and  years,  and  was  greater  in  the  Chesapeake  Bay  Lowlands  than  in  the  Great  Basin. 
Geographic  variation  in  the  extent  to  which  assumptions  of  occupancy  theory  are  met,  and  the 
manner  in  which  such  variation  may  be  associated  with  regional  differences  in  the  biology  of  the 
focal  It  is  possible  that  the  more-heterogeneous  topography  of  the  Great  Basin  mountain  ranges 
in  which  we  worked  leads  to  greater  variation  in  emergence  dates  and  to  more-extensive 
movements  of  individuals  (e.g.,  because  air  flow  along  the  elevational  gradient  can  be  strong) 
than  in  the  relatively  flat  Chesapeake  Bay  Lowlands.  Additionally,  85%  of  the  species  we 
modeled  in  the  Chesapeake  Bay  Lowlands  always  or  sometimes  have  more  than  one  brood  per 
year  (thus  would  be  available  for  sampling  more  consistently,  and  may  be  more  amendable  to 
occupancy  modeling  than  species  with  one  brood  per  year),  compared  to  20%  and  33%  of  the 
species  in  the  central  and  western  Great  Basin,  respectively. 

In  each  ecosystem,  associations  between  one  or  more  covariates  and  py  or  y/i  for  at  least  one- 
third  of  the  species — a  total  of  25  species-by-year  models  across  the  three  ecosystems — either 
could  not  be  estimated  or  seemed  implausibly  large  (i.e.,  >  1 10|).  Examination  of  the  raw  data 
allowed  us  to  identify  potential  causes  of  the  estimation  problems  for  about  2/3  of  these  models. 
In  most  cases,  problems  appeared  to  stem  from  clustering  of  detections  at  one  end  of  the  gradient 
of  values  of  a  covariate.  For  example,  precipitation  values  for  the  transects  in  the  western  Great 
Basin  on  which  S.  zerene  was  detected  in  2014  generally  were  relatively  high.  In  a  few  cases, 
naive  occupancy  may  have  been  too  low  (e.g.,  0.35)  or  too  high  (e.g.,  0.66)  to  allow  estimation 
of  parameters. 

Chesapeake  Bay  Lowlands.  Maximum  single-year  estimates  of  detection  probabilities  (py) 
ranged  from  0.14  to  0.74.  Nectar  was  associated  with  py  of  nine  species  in  one  of  two  years  and 
three  species  in  both  years.  In  all  but  one  case,  as  the  abundance  of  nectar  increased,  py 
increased.  Occupancy  ( y/i )  ranged  from  0.28  to  0.89.  We  identified  covariates  associated  with  y/i 
of  five  species  in  one  year  and  eight  species  in  two  years.  In  the  latter  cases,  one  or  more  of  the 
same  covariates  was  associated  with  yn  of  five  of  the  species  in  both  years.  The  deciduous 
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proportion  of  the  basal  area  of  trees  was  associated  with  i///  of  nine  species  in  2013  and  six 
species  in  2014.  Occupancy  of  one  species  in  2013  and  five  species  in  2014  was  associated  with 
the  number  of  deciduous  stems.  Structural  heterogeneity  was  associated  with  i///  of  four  species 
in  2013  and  two  species  in  2014.  Edge  length  was  associated  with  y/i  °f  one  species  each  in  2013 
and  2014.  In  2013,  the  interaction  between  edge  length  and  structural  heterogeneity  was 
associated  with  y/i  of  one  species. 

In  part,  the  association  between  occupancy  and  either  structural  heterogeneity  or  the  dominance 
or  abundance  of  deciduous  woody  species  may  reflect  that  the  larval  host  plants  of  five  of  the 
species  we  modeled  are  deciduous.  Additionally,  structural  heterogeneity  may  be  a  surrogate 
measure  of  intensity  of  browsing  by  white-tailed  deer.  Areas  in  which  white-tailed  deer  are 
abundant  tend  to  have  relatively  little  understory,  thus  relatively  few  grasses  or  forbs  that  serve 
as  larval  host  plants  or  nectar  sources.  Removal  of  the  understory  also  changes  microclimate  and 
may  expose  immature  life  stages  to  higher  probabilities  of  predation. 

Central  Great  Basin.  Maximum  single-year  py  ranged  from  0.28  to  0.97.  Abundance  of  nectar 
was  associated  with  py  of  three  species  in  2013  only,  two  species  in  2014  only,  and  two  species 
in  both  years.  Abundance  of  mud  was  associated  with  py  of  two  species  in  2013  only,  one 
species  in  2014  only,  and  six  species  in  both  years.  Detection  probability  increased  as  abundance 
of  nectar  or  mud  increased  in  all  but  one  case.  Single-year  y/i  ranged  from  0.44  to  0.98. 
Regression  coefficients  could  not  be  estimated  or  were  >  1 10|  for  one  species  in  1995  and  two 
species  each  in  2013  and  2014.  In  all  other  cases,  one  or  more  covariates  were  associated  with  y/i 
of  all  species  in  each  year.  Occupancy  of  two  species  was  associated  with  elevation  only.  At  least 
some  of  the  covariates  that  were  associated  with  y/i  of  the  other  species  varied  among  years. 
Elevation  (whether  as  a  linear  or  a  quadratic  function)  was  associated  with  y/i  of  six  species  in 
1995,  11  species  in  2013,  and  13  species  in  2014.  Occupancy  of  six  species  in  2013  and  five 
species  in  2014  was  associated  with  mud  abundance.  Precipitation  was  associated  with  y/i  of 
seven  species  in  1995  and  three  species  in  each  of  2013  and  2014.  Nectar  abundance  was 
associated  with  y/i  of  four  species  in  2013  and  five  in  2014.  Occupancy  of  three  species  in  1995, 
four  in  2013,  and  three  in  2014  was  associated  with  terrain  ruggedness. 

Western  Great  Basin.  Maximum  single-year  py  for  a  given  species  ranged  from  0.34  to  0.99. 
Abundance  of  nectar  was  associated  with  py  of  two  species  in  2014  only  and  four  species  in  both 

2013  and  2014.  Abundance  of  mud  was  associated  with  py  of  four  species  in  2013  only,  two  in 

2014  only,  and  three  in  both  years.  Single-year  y/i  ranged  from  0.31  to  0.78.  Regression 
coefficients  and  effect  sizes  could  not  be  estimated  or  were  >  |10|  for  one  species  each  in  2013 
and  2014.  In  all  other  cases,  one  or  more  covariates  were  associated  with  y/i  of  all  species  in  both 
years.  A  single  covariate — elevation  or  precipitation — was  associated  with  y/i  of  three  species  in 
both  2013  and  2014.  At  least  some  of  the  covariates  that  were  associated  with  y/i  of  the  other 
species  varied  between  years.  Elevation  (whether  as  a  linear  or  a  quadratic  function)  was 
associated  with  y/i  of  nine  species  in  2013  and  six  species  in  2014.  Precipitation  was  associated 
with  y/i  of  eight  species  in  2013  and  seven  in  2014.  Occupancy  of  five  species  in  either  year  was 
associated  with  terrain  ruggedness.  Mud  abundance  was  associated  with  y/i  of  one  species  in 
2013  and  one  species  in  both  years.  Nectar  abundance  was  associated  with  y/i  of  two  species  in 
2014. 
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Either  elevation  or  precipitation  were  associated  with  occupancy  of  all  species  in  the  central  and 
western  Great  Basin.  Strong  relations  between  elevation  and  occupancy  were  consistent  with 
previous  work  in  these  ecosystems  (e.g.,  Fleishman  et  al.  1998,  1999,  2001a,  b).  In  the  central 
Great  Basin,  about  twice  as  many  associations  between  occupancy  and  precipitation  were 
observed  in  1995,  a  water  year  with  unusually  high  precipitation,  than  in  2013  or  2014.  The 
direction  of  the  associations  between  occupancy  and  precipitation  was  inconsistent  among 
species.  Although  extreme  winters  can  increase  mortality  of  overwintering  life  stages,  especially 
larvae  (Douglas  1986,  Dennis  1993),  precipitation  may  forestall  senescence  of  host  plants  and 
nectar  sources. 

Our  estimates  of  abundance  of  nectar  and  mud  were  coarse;  they  were  intended  to  be  rapid  and 
fairly  repeatable  among  observers.  Nevertheless,;?//  in  one  or  more  years  was  associated  with 
abundance  of  nectar  or  mud  for  92%  of  the  species  in  the  Chesapeake  Bay  Lowlands,  all  of  the 
species  in  the  central  Great  Basin,  and  93%  of  the  species  in  the  western  Great  Basin.  In  almost 
all  cases,  py  increased  as  the  abundance  of  these  resources  increased.  These  relations  are 
consistent  with  field  experience.  For  example,  species  in  the  Papilionidae  and  Polyommatinae 
often  are  detected  at  mud.  However,  to  the  best  of  our  knowledge,  these  relations  have  not 
previously  been  quantified  at  the  assemblage  level.  Detection  probabilities  tended  to  be  relatively 
high  for  species  that  are  abundant,  have  limited  vagility  and  circumscribed  habitat,  or  are 
conspicuous.  Detection  probability  tended  to  be  associated  more  consistently  than  occupancy 
with  abundance  of  nectar  or  mud.  Continuous  estimates  of  sugar  mass  may  be  more  closely 
associated  with  occupancy  than  categorical  abundance  of  nectar  (Pavlik  et  al.  in  review). 

We  modeled  seven  species  that  occurred  in  both  the  central  and  western  Great  Basin,  but 
covariates  associated  with  py  and  y/i  of  those  species,  or  the  direction  of  association  with  the 
same  covariate,  often  differed  among  ecosystems.  This  variation  may  reflect  that  our  study 
locations  in  the  central  and  western  Great  Basin  are  within  different  centers  of  differentiation  or 
zoogeographic  regions  (Austin  and  Murphy  1987).  In  many  cases,  different  subspecies  occupy 
the  two  regions,  and  their  local  ecology  may  be  sufficiently  distinct  to  affect  covariate 
associations. 

Occupancy  rarely  has  been  estimated  for  a  high  proportion  of  the  species  within  an  assemblage 
of  butterflies,  and  our  work  highlighted  some  of  the  applications  in  which  occupancy  estimation 
has  limitations.  Single-species  occupancy  models  often  are  applied  to  many  species  that  were 
sampled  simultaneously.  These  situations  may  require  or  lead  to  the  assumption  that  sampling  of 
each  modeled  species  was  sufficiently  and  equally  robust.  But  given  typical  financial  and 
logistical  constraints,  it  rarely  is  tractable  to  sample  an  assemblage,  especially  one  that  is  highly 
dynamic,  with  a  design  that  is  ideal  for  estimating  occupancy  of  the  majority  of  individual 
species.  For  example,  assemblage-level  surveys  generally  encompass  areas  that  are  unlikely  to  be 
occupied  by  a  given  species,  which  may  complicate  parameter  estimation.  Moreover,  if 
organisms  are  sufficiently  mobile  that  the  occupancy  status  of  fairly  small  sites  is  likely  to 
change  between  surveys,  as  is  the  case  with  adult  butterflies,  inferences  about  the  environmental 
variables  that  are  associated  with  py  and  y/i  may  become  biased  (e.g.,  Hayes  and  Monfils  2015). 
Our  work  both  elucidates  trade-offs  among  application  of  occupancy  models  to  multiple  co¬ 
occurring  species  of  butterflies  and  highlights  a  number  of  novel  ecological  relations,  especially 
the  extent  to  which  detection  probabilities  may  relate  to  ephemeral  resources. 
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Use  of  occupancy  as  a  surrogate  measure  of  species  richness:  environmental  associations 
with  multiple  states  of  anuran  occupancy 

Results  of  this  ongoing  analysis  improve  understanding  of  assessment  and  monitoring  of  anurans 
at  multiple  spatial  and  temporal  scales  and  provide  information  for  projecting  responses  of 
anurans  to  environmental  change.  We  first  applied  multiple- state  models  to  anurans  because  the 
methods  were  highly  consistent  with  the  manner  in  which  anuran  data  were  collected — a 
categorical  calling  index  rather  than  an  estimated  count  (number  of  animals).  Given  that  our 
application  to  date  generally  has  met  technical  criteria  for  success  (e.g.,  most  models  converged 
and  variances  were  not  extreme)  and  yielded  useful  biological  inferences,  we  believe  there  is 
considerable  potential  to  develop  abundance  classes  for  birds  and  butterflies  and  apply  the 
models  to  additional  data  collected  as  part  of  the  current  project. 

We  detected  Fowler’s  toads,  Cope’s  gray  tree  frogs,  and  spring  peepers  at  64,  73,  and  71  of  108 
wetlands,  respectively.  The  calling  index  of  Fowler’s  toads,  Cope’s  gray  tree  frogs,  and  spring 
peepers  was  high  during  >  1  survey  at  26,  41,  and  49  wetlands,  respectively. 

Estimates  of  detection  probability  for  anurans  at  low  abundance  in  a  single  year  ranged  from 
0.01  ±  0.04  (SE)  (Fowler’s  toad,  2015)  to  0.90  ±  0.05  (spring  peeper,  2013)  (Table  5).  Estimates 
of  detection  probability  for  anurans  at  high  abundance  in  a  single  year  ranged  from  0.23  ±  0.04 
(Fowler’s  toad,  2015)  to  1.0  (spring  peeper,  2012  and  2013).  Multiple-year  estimates  of 
occupancy  at  low  abundance  fell  between  0.51  ±  0.07  (spring  peeper)  and  0.80  ±  0.06  (Fowler’s 
toad),  whereas  estimates  of  occupancy  at  high  abundance  fell  between  0.40  ±  0.10  (Fowler’s 
toad)  and  0.92  ±  0.05  (spring  peeper).  The  probability  of  any  of  the  species  colonizing  a 
previously  vacant  wetland,  whether  in  a  low  or  a  high  abundance  state,  did  not  exceed  0.12 
(Table  5),  whereas  the  probability  of  extirpation,  remaining  at  low  abundance,  or  transitioning 
from  high  to  low  abundance  was  a  minimum  of  0.80.  The  probability  of  remaining  at  high 
abundance  fell  between  0.18  ±  0.09  (Cope’s  grey  tree  frog)  and  0.48  ±  0.10  (Fowler’s  toad). 

Although  scale  is  referenced  often  in  ecology,  models  that  include  covariates  measured  at 
multiple  extents  are  relatively  uncommon.  Occupancy  of  pond-breeding  anurans  was  associated 
with  variables  at  multiple  spatial  extents  that  are  relevant  to  different  stages  of  their  life  history. 
The  association  between  occupancy  and  environmental  variables  at  the  smallest  extent,  which 
corresponded  to  breeding,  was  poor  to  fair  (AUC  0.60-0.75),  and  was  equal  or  weaker  than 
associations  between  occupancy  and  environmental  variables  at  extents  that  corresponded  to 
migration  (AUC  0.70-0.84)  and  dispersal  (AUC  0.64-0.74).  Models  that  included  covariates  at 
multiple  extents  generally  explained  more  variation  in  occupancy  than  models  that  included 
covariates  at  one  extent  (AUC  0.71-0.92). 

Both  single-state  and  multiple-state  models  of  occupancy  of  Fowler’s  toads  suggested  that 
occupancy  was  strongly  associated  with  the  magnitude  of  forest  aggregation  within  400  m  of  the 
wetland  (positive),  the  percentage  of  upland  forest  cover  within  400  m  (positive),  the  density  of 
highways  with  high  traffic  volume  within  4  km  (negative),  effective  mesh  size  within  1.5  km 
(positive),  the  percentage  of  wetland  cover  within  3.5  km  (positive),  and  the  percentage  of 
impervious  surface  within  2.5  km  (negative).  Percentage  of  canopy  may  serve  as  a  proxy 
measure  of  shoreline  vegetation  structure. 
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Table  5.  Preliminary  results  of  multiple- season,  multiple-state  analyses  of  the  occupancy  of  anurans  in  the  Chesapeake  Bay  Lowlands. 
Psi,  occupancy.  State  1,  low  abundance;  state  2=  high  abundance.  SE,  standard  error. _ _ 


Spring  peeper 

Cope’s  grey  tree  frog 

Fowler’s  toad 

Estimate 

Estimate 

SE 

Estimate 

SE 

Estimate 

SE 

Occupancy,  low  abundance 

0.51 

0.07 

0.64 

0.07 

0.80 

0.06 

Transition  from  absent  to  low  abundance 

0.01 

0.01 

0.12 

0.05 

0 

NA 

Remain  at  low  abundance 

0.92 

0.06 

0.90 

0.07 

1.0 

NA 

Transition  from  low  abundance  to  absent 

1.0 

NA 

1.0 

NA 

0.94 

0.04 

Occupancy,  high  abundance 

0.92 

0.05 

0.75 

0.11 

0.40 

0.10 

Transition  from  absent  to  high  abundance 

0 

0 

0.09 

0.20 

0 

0 

Remain  at  high  abundance 

0.45 

0.11 

0.18 

0.09 

0.48 

0.10 

Transition  from  high  to  low  abundance 

0.82 

0.64 

0.80 

0.07 

0.93 

0.06 

Detection  probability,  low  abundance,  2011 

0.49 

0.21 

0.35 

0.09 

0.06 

0.02 

Detection  probability,  low  abundance,  2012 

0.27 

0.16 

0.13 

0.05 

0.11 

0.04 

Detection  probability,  low  abundance,  2013 

0.90 

0.05 

0.17 

0.06 

NA 

NA 

Detection  probability,  low  abundance,  2014 

0.66 

0.13 

0.10 

0.03 

0.08 

0.05 

Detection  probability,  low  abundance,  2015 

0.48 

0.15 

0.24 

0.06 

0.01 

0.04 

Detection  probability,  high  abundance,  2011 

0.97 

0.02 

0.47 

0.04 

0.39 

0.07 

Detection  probability,  high  abundance,  2012 

1.0 

0 

0.41 

0.06 

0.30 

0.04 

Detection  probability,  high  abundance,  2013 

1.0 

0 

0.58 

0.04 

0.38 

0.04 

Detection  probability,  high  abundance,  2014 

0.96 

0.02 

0.47 

0.08 

0.34 

0.04 

Detection  probability,  high  abundance,  2015 

0.97 

0.02 

0.54 

0.05 

0.23 

0.04 
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Both  single-state  and  multiple- state  models  of  occupancy  of  Cope’s  gray  tree  frogs  suggested 
that  occupancy  was  strongly  associated  with  whether  the  wetland  was  open  water  or  had 
emergent  vegetation,  the  percentage  of  upland  forest  cover  within  1  km  (positive),  the  magnitude 
of  forest  aggregation  within  1  km  (positive),  the  percentage  of  impervious  surface  within  3.5  km, 
the  density  of  highways  with  high  traffic  volume  within  500  m  (negative),  effective  mesh  size 
within  5  km  (positive),  and  the  percentage  of  wetland  cover  within  5  km.  The  species  was  more 
strongly  associated  with  wetlands  than  with  ponds. 

Both  single-state  and  multiple- state  models  of  occupancy  of  spring  peepers  suggested  that 
occupancy  was  strongly  associated  with  the  percentage  of  canopy  cover  around  a  wetland 
(positive),  the  magnitude  of  forest  aggregation  and  the  percentage  of  forest  cover  within  600  m 
(positive),  the  percentage  of  impervious  surface  within  4.5  km  (negative),  effective  mesh  size 
within  1.5  km  or  4.5  km  (positive),  and  wetland  density  within  3.5  km  (positive).  The  single¬ 
state  analysis  suggested  that  occupancy  was  positively  associated  with  the  length  of  the 
perimeter  of  the  wetland,  whereas  the  multiple- state  analysis  suggested  that  occupancy  was 
associated  with  wetland  area  (positive)  and  highway  density  within  5  km  (negative). 

Many  anuran  populations  are  believed  to  function  as  metapopulations  (Smith  and  Green  2005). 
Although  the  process  of  anuran  dispersal  is  poorly  understood  (Semlitsch  2008),  occupancy  is 
negatively  related  to  urbanization  (Rubbo  and  Kiesecker  2005),  roads  (Eigenbrod  et  al.  2008), 
and  distance  to  neighboring  wetlands  (Scherer  et  al.  2012).  The  distance  at  which  many  of  the 
migration-level  and  dispersal-level  covariates  were  associated  with  occupancy  was  the  maximum 
quantified,  suggesting  that  dispersal  may  be  occurring  across  areas  larger  than  we  defined. 

Our  preliminary  results  suggest  that  populations  of  Fowler’s  toads,  Cope’s  gray  tree  frogs,  and 
spring  peepers  in  the  Chesapeake  Bay  Lowlands  generally  are  stable,  but  if  these  anurans  are 
extirpated  from  wetlands  in  the  region,  the  wetlands  are  unlikely  to  be  recolonized. 


Responses  of  songbird  occupancy  to  environmental  change 

Results  of  this  analysis  provide  information  that  facilitates  projection  of  the  effects  of  land  use 
and  management  actions  on  native  birds  in  the  Chesapeake  Bay  Lowlands.  Although 
urbanization  is  likely  to  continue  throughout  the  ecosystem,  densities  of  white-tailed  deer  can  be 
managed  via  culling. 

To  the  best  of  our  knowledge,  we  conducted  the  first  statistically  rigorous,  multiple-year 
investigation  of  regional  relations  between  intensity  of  deer  browsing  and  densities  of  individual 
species  and  guilds  of  birds.  Our  results  were  consistent  with  a  growing  body  of  evidence  that 
deer  affect  the  quality  of  habitat  for  birds,  especially  woodland  songbirds  that  depend  on  foliage 
below  the  browse  line  (e.g.,  McShea  and  Rappole  2000,  Cardinal  et  al.  2012).  Areas  in  which 
densities  of  deer  are  high  have  little  vegetation  cover  below  the  browse  line  (Martin  et  al.  2010), 
and  direct  depredation  of  songbird  nests  by  deer  has  been  documented  (Luepold  et  al.  2015). 

We  estimated  detection  probability  as  a  nuisance  parameter  to  adjust  pellet  and  bird  densities. 
The  probability  of  detecting  deer  pellets  was  significantly  lower  in  coastal  Virginia  than  inland 
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Virginia.  All  detection  probabilities  of  individual  species  and  guilds  of  birds  were  <  0.34. 
Detection  probabilities  at  the  guild  level  and  at  the  level  of  species  within  that  guild  did  not 
differ  in  coastal  Virginia.  The  probability  of  detecting  species  that  nest  in  understory  shrubs  and 
glean  understory  foliage  did  not  differ  between  regions.  However,  the  probability  of  detecting 
the  other  two  guilds  was  lower  in  inland  than  coastal  Virginia. 

Relatively  high  intensity  of  deer  use  in  coastal  Virginia  was  correlated  with  substantial  levels  of 
local  forest  fragmentation,  which  is  associated  with  high-quality  deer  habitat.  The  mean  extent  of 
forest-road,  forest-grassland,  forest-row  crop,  and  forest-rural  edges  was  232%,  346%,  176%, 
and  252%  higher,  respectively,  in  coastal  Virginia  than  inland  Virginia.  Deer  use  in  2012  was 
greater  in  coastal  Virginia  (2014  pellets  ha'1,  range  0-28,193  ha"1)  than  inland  Virginia  (0  pellets 
ha'1,  range  0-19,600  pellets  ha'1)  (W  =  5,462,  P  =  0.01).  We  detected  no  pellets  in  38%  of  the 
sites  in  coastal  Virginia  and  57%  of  the  sites  in  inland  Virginia.  The  mean  number  of  deer 
harvested  in  coastal  Virginia  from  2004-2013  (4.73,  Cl  4.45-5.00)  was  greater  than  in  inland 
Virginia  (2.99,  Cl  2.79-3.19)  (ft  =  11.1,  P  <  0.001). 

Deer  use  and  the  density  of  shrub-nesting  foliage  gleaners  in  coastal  Virginia  was  negatively 
correlated  (rs  =  -0.35,  P  <  0.01)  (Figure  10).  In  coastal  Virginia,  deer  use  and  the  density  of 
canopy-nesting  or  cavity-nesting  species  that  feed  in  the  canopy  or  sally  aerially  was  positively 
correlated  (rs  =  0.27,  P  =  0.03).  Deer  use  and  the  density  of  guilds  in  inland  Virginia  was  not 
significantly  correlated.  In  coastal  Virginia,  deer  use  was  negatively  correlated  with  the  density 
of  White-eyed  Vireos  ( Vireo  griseus).  Red-eyed  Vireos  ( Vireo  olivaceus ),  Black-and-white 
Warblers  ( Mniotilta  varia ),  Hooded  Warblers  ( Setophaga  citrina).  Prairie  Warblers  ( Setophaga 
discolor ),  and  Scarlet  Tanagers  ( Piranga  olivacea ).  Prairie  Warblers,  which  breed  in  forests  with 
open  canopy  and  shrubby  undergrowth  (Nolan  et  al.  1999),  are  classified  as  a  species  of  highest 
conservation  concern  at  the  continental  level  by  Partners  in  Flight  (Rosenberg  et  al.  2016).  Local 
populations  of  Hooded  Warblers,  which  nest  in  understory  shrubs,  have  declined  considerably  as 
the  shrub  layer  has  disappeared  (Chiver  et  al.  2011). 

Deer  use  was  positively  correlated  with  the  density  of  Eastern  Wood-Pewees  ( Contopus  virens), 
Acadian  Flycatchers  ( Empidonax  virescens ),  Blue  Jays  ( Cyanocitta  cristata),  Carolina 
Chickadees  ( Poecile  carolinensis),  and  White -breasted  Nuthatches  (, Sitta  carolinensis).  Some  of 
these  correlations  may  reflect  overlap  in  the  distribution  of  habitat  quality  for  deer  and  some 
species  of  birds.  For  example,  both  Carolina  Chickadees  and  deer  thrive  in  wooded  suburbs 
(Mostrom  et  al.  2002). 
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Figure  10.  Relations  between  density  of  birds  in  different  guilds  and  relative  intensity  of  deer 
use  (density  of  deer  pellets).  Left,  shrub-nesting  foliage  gleaners;  middle,  species  that  nest  in  the 
canopy  and  feed  on  open  ground;  right,  middle,  canopy-nesting  or  cavity-nesting  species  that 
feed  in  the  canopy  or  sally  aerially. 


Partitioning  drivers  of  beta  diversity 

Results  of  this  analysis  may  increase  the  ability  to  develop  effective  and  cost-efficient  methods 
for  assessing  and  monitoring  native  birds  and  butterflies  in  multiple  subregions  of  the  Great 
Basin  across  space  and  time.  Assessment  and  monitoring  of  native  faunas  is  likely  to  become 
increasingly  relevant  if  removal  of  conifers,  intended  to  increase  habitat  quality  for  Greater  Sage- 
Grouse,  is  implemented  as  planned. 

Temporal  patterns.  Our  preliminary  results  indicated  that  bird  species’  identities  were  more 
consistent  through  time  than  would  be  expected  under  a  null  model,  which  suggests  that 
distributions,  although  variable,  are  associated  with  abiotic  or  biotic  environmental  attributes. 
More  variation  in  beta  diversity  was  attributable  to  turnover  than  to  nestedness — increases  in 
species  richness  do  not  simply  represent  addition  of  species  to  a  core  group.  However,  observed 
nestedness  was  generally  higher  than  under  a  null  model,  which  is  consistent  with  the  fact  that 
the  species  pool  is  relatively  similar  across  the  region. 

Temporal  patterns  in  turnover  and  nestedness  of  butterfly  assemblages  were  similar  to  those  of 
bird  assemblages:  turnover  was  lower  than  expected  and  nestedness  higher  than  expected  under  a 
null  model.  Levels  of  nestedness  of  butterflies  were  equal  to  or  greater  than  those  of  birds, 
whereas  turnover  was  lower. 

Little  of  the  variation  in  temporal  turnover  and  nestedness  was  explained  by  zoogeographic 
region  or  spatial  resolution,  although  several  weak  associations  were  apparent.  Temporal 
turnover  and  nestedness  of  bird  assemblages  were  higher  in  the  central  Great  Basin  than  in  the 
western  Great  Basin,  and  at  the  level  of  points  than  at  the  level  of  canyons.  Turnover  of  butterfly 
assemblages  was  lower  in  the  central  Great  Basin  than  in  the  western  Great  Basin,  but  did  not 
differ  among  spatial  resolutions.  The  level  of  nestedness  of  butterfly  assemblages  did  not  differ 
between  regions  or  spatial  resolutions. 

Spatial  patterns.  Bird  species’  identities  again  were  more  consistent  in  space  than  would  be 
expected  under  a  null  model,  and  were  more  consistent  among  canyons  within  mountain  ranges 
than  among  points  within  canyons.  Spatial  nestedness  in  bird  assemblages  was  often  higher  than 
expected  under  a  null  model,  and  was  lower  in  the  western  than  in  the  central  Great  Basin.  The 
majority  of  spatial  variation  in  bird  assemblages  was  due  to  turnover  rather  than  nestedness. 

Spatial  turnover  in  butterfly  assemblages  was  often  lower  than  expected  under  a  null  model,  and 
was  higher  in  the  western  than  in  the  central  Great  Basin.  This  may  reflect  that  there  is  more 
biogeographic  differentiation  in  our  western  Great  Basin  study  areas  than  in  our  study  areas  in 
the  central  Great  Basin.  The  avifauna  of  the  east  slope  of  the  Sierra  Nevada  includes  several 
species  that  typically  do  not  occur  further  east  in  the  Great  Basin,  and  vice  versa.  Spatial 
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turnover  was  lower  in  butterfly  assemblages  than  in  bird  assemblages.  Spatial  nestedness  in 
butterfly  assemblages  often  was  higher  than  in  bird  assemblages.  The  latter  is  consistent  with  the 
fact  that  although  butterflies  have  fairly  specific  habitat  requirements  as  larvae,  their  adult  habitat 
requirements  may  be  more  similar  than  that  of  birds. 
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Conclusions  and  Future  Research  and  Implementation 


We  achieved  the  objectives  of  the  statement  of  need  to  which  the  project  responded  by  providing 
knowledge  that  improves  understanding  of  assessment,  monitoring,  and  management  of  diverse 
native  species  at  multiple  spatial  and  temporal  scales.  This  knowledge  and  understanding  can 
directly  inform  development  of  effective  and  cost-efficient  methods  for  assessment,  monitoring, 
and  projection  of  the  effects  of  land  use  and  management  actions.  We  empirically  tested  diverse 
aspects  of  theory  related  to  assessment  and  monitoring  of  native  species,  including  the  fact  that 
habitat  is  species-specific,  and  therefore  different  species  will  have  different  responses  to 
environmental  change.  Our  work  provides  information  on  the  extent  to  which  theory  is 
applicable  to  multiple  taxonomic  groups  that  often  are  the  focus  of  assessment  and  monitoring, 
and  how  such  applicability  varies  within  extensive  ecosystems  and  among  ecosystems. 

In  some  cases,  we  met  the  project  objectives  in  ways  that  were  different  from  what  we  initially 
expected.  For  example,  we  anticipated  that  we  would  incorporate  estimates  of  occupancy  into 
our  analyses  of  the  explanatory  and  predictive  ability  of  indicator-species  models.  However,  we 
discovered  through  analyses  of  birds  and  butterflies  that  because  detection  probabilities  for  many 
species  were  below  0.3,  we  often  could  not  estimate  occupancy  for  more  than  50%  of  the  species 
detected  in  each  ecosystem  and  year.  As  a  consequence,  our  response  variable  in  hierarchical 
modeling  of  species  richness,  indicator  species,  and  beta  diversity  was  naive  occupancy  rather 
than  detection-weighted  occupancy.  These  results  provide  empirical  evidence  that  some 
applications  of  occupancy  models  that  have  been  suggested  in  the  literature  may  not  be  effective 
in  practice,  especially  for  organisms  with  relatively  dynamic  short-term  distributions  or  in 
ecosystems  in  which  weather  and  climate  are  relatively  variable  and  unpredictable. 


Extents  of  sampling  and  management 

Hierarchical  estimates  of  species  richness  and  analyses  of  beta  diversity  highlighted  differences 
in  the  spatial  extent  at  which  birds  and  butterflies  appear  to  respond  to  their  environment.  In  the 
mountains  of  the  Great  Basin,  canyons  appear  to  be  a  more  biogeographically  meaningful  spatial 
unit  for  birds  than  relatively  small  points  that  typically  are  used  as  sampling  units  for  that 
taxonomic  group.  Therefore,  it  may  be  worthwhile  to  sample  birds  in  multiple  locations  in  a 
given  canyon,  and  to  include  canyon  as  a  covariate  in  analyses,  rather  than  treating  points  as 
replicates  independent  of  the  canyon  within  which  they  are  embedded. 

In  our  hierarchical  analyses  of  species  richness,  the  proportion  of  variance  in  canyon-level 
cumulative  species  richness  of  birds  that  we  were  able  to  explain  as  a  function  of  environmental 
covariates  was  at  least  double  the  proportion  we  could  explain  at  the  point  level  in  both  the 
central  and  western  Great  Basin.  Moreover,  a  combination  of  canyon-level  and  point-level 
processes  explained  more  variation  in  cumulative  species  richness  of  birds  than  did  point-level 
processes  alone.  These  results  suggest  that  the  effects  of  land  use  or  management  in  one  area  of  a 
canyon,  especially  if  related  to  riparian  areas,  may  have  indirect  effects  on  birds  that  breed 
throughout  the  canyon. 


70 


By  contrast,  hierarchical  analyses  of  species  richness  allowed  us  to  explain  more  than  0.70  of  the 
transect-level  variation  in  cumulative  species  richness  of  butterflies  in  both  the  central  and 
western  Great  Basin.  Similarly,  in  analyses  of  beta  diversity  in  the  central  and  western  Great 
Basin,  temporal  turnover  and  nestedness  of  bird  assemblages  was  greater  at  the  level  of  points 
than  the  level  of  canyons,  whereas  temporal  beta  diversity  of  butterflies  did  not  differ  between 
spatial  resolutions.  These  results  suggest  that  management  of  butterfly  habitat  may  be  effective  at 
relatively  small  extents,  whereas  management  of  bird  habitat  may  need  to  be  undertaken  at  larger 
extents. 

The  extent  to  which  processes  at  multiple  extents  are  associated  with  local  occupancy  also  was 
apparent  in  preliminary  results  from  multiple- state  models  of  anuran  occupancy.  Although  many 
populations  appeared  to  be  stable,  we  found  little  evidence  of  colonization  by  anurans,  which 
suggests  that  if  anurans  are  extirpated  from  a  given  natural  or  anthropogenic  wetland  in  the 
Chesapeake  Bay  Lowlands,  that  location  is  unlikely  to  be  recolonized. 


Estimation  of  occupancy  versus  species  richness 

We  noted  above  that  although  the  difference  in  estimates  of  bird  occupancy  that  were  based  on 
5-min  and  8-min  counts  in  the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin 
was  slight,  it  might  be  sufficient  to  affect  estimates  of  species  richness  of  birds  that  are  based  on 
occupancy  models  (Iknayan  et  al.  2014),  especially  because  sample  sizes  were  too  small  to 
estimate  occupancy  for  a  substantial  proportion  of  the  avifauna.  We  conducted  a  simple  test  of 
the  effects  of  count  duration  on  naive  estimates  of  species  richness  of  birds  (i.e.,  simple  counts 
not  standardized  by  area  to  which  sampling  was  extrapolated,  number  of  individuals  sampled,  or 
detection  probability)  at  nested  spatial  extents.  This  relation  has  been  examined  in  many  studies, 
but  evidence  is  equivocal  and  likely  reflects  variation  among  land-cover  types  (Petit  et  al.  1995, 
Shiu  and  Lee  2003),  ecosystems,  sampling  designs,  and  analysis  methods.  The  number  of  species 
detected  increases  as  count  duration  increases,  sometimes  substantially  even  at  1  min  increments 
(Thompson  and  Schwalbach  1995),  but  perhaps  less  so  when  estimates  of  species  richness  are 
based  on  multiple  counts  throughout  the  season  rather  than  one  count  (Buskirk  and  McDonald 
1995).  In  general,  it  is  not  uncommon  for  ca.  65-90%  of  detections  of  a  given  species  (Reidy  et 
al.  2011),  of  all  species  (Tegler-Amones  et  al.  2012),  or  of  the  estimated  mean  number  of  species 
to  be  recorded  within  the  first  5  min  of  a  count  (Esquivel  and  Peris  2008).  By  contrast,  the 
percentage  of  mean  deviance  explained  by  generalized  additive  models  of  distributions  of  bird 
species  in  southwestern  France  that  were  based  on  either  5  min  or  10  min  counts  was  the  same 
(25%);  increasing  the  duration  to  20  min  increased  the  mean  deviance  explained  to  28% 
(Bonthoux  and  Balent  2012).  We  aimed  not  only  to  better  understand  sampling  trade-offs  in  our 
own  study  systems  but  to  contribute  to  ongoing,  rigorous  assessments  of  field  methods,  including 
methods  associated  with  counts  of  organisms. 

We  calculated  species  richness  (number  of  species  counted,  not  weighted  by  detection 
probability)  of  birds  in  2012  and  2013  at  the  level  of  points.  In  the  western  and  central  Great 
Basin  we  also  calculated  species  richness  at  two  nested  levels:  canyons  or  polygons  defined  by 
previous  or  planned  management  treatments,  and  mountain  ranges.  Where  species  richness 
estimated  on  the  basis  of  5  min  and  8  min  counts  differed,  we  identified  the  additional  species 
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that  were  detected  from  5-8  min.  We  compared  mean  species  richness  at  the  level  of  points  in 
the  Chesapeake  Bay  Lowlands  and  central  and  western  Great  Basin,  and  at  the  level  of  canyons 
in  the  central  and  western  Great  Basin,  with  one-tailed  paired  /-tests  in  StatPlus  (v.5.8, 
AnalystSoft,  Alexandria,  Virginia).  We  used  counts  for  this  comparison  because  the  proportion 
of  species  for  which  sample  sizes  were  insufficient  to  model  occupancy  was  high. 

Data  were  sufficient  to  estimate  occupancy  for  26  of  the  62  species  in  the  Chesapeake  Bay 
Lowlands.  Data  were  sufficient  to  estimate  occupancy  for  18  of  the  79  species  detected  in  the 
central  Great  Basin  in  2012  (23%)  and  26  of  the  80  species  detected  in  2013  (33%).  We 
estimated  occupancy  for  25  of  the  87  species  detected  in  the  western  Great  Basin  in  2012  (29%) 
and  29  of  the  86  species  detected  in  2013  (34%). 


In  all  three  ecosystems,  both  point-level  and  canyon-level  estimates  of  species  richness  based  on 
8  min  counts  consistently  were  higher  than  those  based  on  5  min  counts,  and  all  differences  were 
statistically  significant  (Table  4).  There  was  no  clear  pattern  in  the  identity  or  attributes  of 
species  that  were  detected  from  5-8  min  in  a  given  location  or  year. 


Table  4.  Naive  estimates  of  species  richness  of  birds  at  the  point  level  and,  in  the  central  and 
western  Great  Basin,  at  the  canyon  level  on  the  basis  of  points  sampled  for  5  and  8  min.  All  tests 
statistically  significant  (p  <  0.001).  SD,  standard  deviation. 


mean  SD  df  t 

Point  level 


Chesapeake  Bay  Lowlands 


2012  5  min 

13.82 

3.10 

130 

18.65 

8  min 

15.92 

3.31 

2013  5  min 

13.68 

3.14 

130 

15.84 

8  min 

15.61 

3.47 

central  Great  Basin 

2012  5  min 

7.08 

3.18 

295 

16.89 

8  min 

8.23 

3.63 

2013  5  min 

7.88 

3.24 

313 

17.82 

8  min 

9.08 

3.71 

western  Great  Basin 

2012  5  min 

10.08 

3.05 

157 

15.30 

8  min 

11.64 

3.48 

2013  5  min 

10.44 

3.11 

157 

15.35 

8  min 

11.89 

3.35 

Canyon  level 
central  Great  Basin 

2012  5  min  24.59  7.68  26  6.91 
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8  min 

26.81 

8.58 

2013  5  min 

25.28 

8.04 

28 

9.44 

8  min 

27.69 

8.29 

western  Great  Basin 

2012  5  min 

31.77 

6.37 

12 

5.59 

8  min 

34.54 

5.92 

2013  5  min 

33.08 

5.53 

12 

5.84 

8  min 

35.00 

5.34 

Inferences  from  occupancy  models 


The  primary  focus  of  our  occupancy  models  was  elucidating  species-environment  relations 
rather  than  precise  estimation  of  occupancy  per  se.  In  that  context,  if  one  must  decide  between 
increasing  the  number  of  sample  units  or  increasing  the  number  of  visits,  we  believe  that 
increasing  the  number  of  sample  units  remains  more  effective.  Additionally,  estimation  of  use 
may  be  equally  or  more  informative  than  estimation  of  occupancy  for  exploring  species- 
environment  relations.  Given  that  few  if  any  of  our  systems  appear  to  have  been  closed,  it  is 
likely  that  our  models  largely  estimated  use. 

Our  mechanistic  research  on  Wood  Thrushes  ( Hylocichla  mustelina )  in  the  Chesapeake  Bay 
Lowlands  (Jirinec  et  al.  2016)  highlighted  the  extent  to  which  the  availability  of  species  for 
sampling  is  likely  to  change  within  a  season,  even  for  species  that  are  territorial.  Jirinec  et  al. 
(2016)  tracked  37  radio-tagged  males  and  used  the  resulting  data  to  construct  95%  kernel  diurnal 
home  ranges.  They  also  tagged  and  tracked  the  female  mates  in  10  home  ranges.  Of  74  male 
night  roosts,  31%  were  located  outside  diurnal  home  ranges. 

Our  work  highlighted  the  extent  to  which  abundance  of  nectar  and  mud  may  be  associated  with 
detection  probabilities  and  occupancy  or  use  of  butterflies.  Nectar  and  mud  generally  have  been 
underemphasized  in  assessment  and  monitoring  programs  for  butterflies,  whereas  more  emphasis 
has  been  placed  on  larval  hostplants.  Previous  laboratory  and  controlled  experiments 
demonstrated  that  some  species  of  butterflies  may  prefer  sucrose  to  other  sugars.  However,  sugar 
use  has  not  been  quantified  for  an  assemblage  of  butterflies  in  the  field.  The  abundance  of 
nectar-producing  plants,  and  the  volume  and  concentration  of  the  nectar  in  those  plants,  may 
peaks  in  the  initial  years  following  a  major  disturbance  such  as  fire,  facilitating  investigation  of 
nectar  use  by  butterflies.  In  2014  and  2015,  we  surveyed  butterflies  and  vegetation  within  the 
boundary  of  the  Rim  Fire  (Stanislaus  National  Forest,  Tuolumne  County,  California),  a  major 
fire  that  occurred  in  2013.  We  quantified  the  masses  of  glucose,  fructose,  and  sucrose  for  all 
plant  species  on  which  we  observed  butterflies  feeding.  We  tested  whether  intensity  of  use  of 
each  20  nectar  sources  (the  number  of  butterflies  observed  taking  nectar  from  each  source  across 
both  years)  was  associated  with  the  total  sugar  mass,  mass  of  sucrose,  or  relative  proportion  of 
sucrose.  We  found  no  evidence  that  intensity  of  butterfly  use  was  associated  with  sugar  mass, 
mass  of  sucrose,  or  the  relative  proportion  of  sucrose  (Pavlik  et  al.  unpublished  manuscript). 
Instead,  butterflies  appeared  to  use  indiscriminately  any  nectar  sources  that  were  available  to 
them.  The  difference  between  apparent  sugar  preferences  in  the  laboratory  and  in  the  field  may 
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be  explained  in  part  by  resource  availability,  and  may  change  as  vegetation  succession 
progresses. 

Our  work  allowed  us  to  draw  inference  to  how  species  richness  and  occupancy  of  birds  may 
respond  to  environmental  changes  in  the  Chesapeake  Bay  Lowlands  and  central  and  western 
Great  Basin  that  largely  can  be  controlled  by  management.  For  example,  we  were  able  to  gauge 
the  extent  to  which  given  levels  of  fragmentation  of  deciduous  forest  in  the  Chesapeake  Bay 
Lowlands,  and  given  densities  of  white-tailed  deer,  may  positively  or  negatively  affect  density  of 
individual  species  and  functional  groups  of  songbirds,  including  species  that  are  declining  across 
the  region.  Ongoing  analyses  of  responses  of  functional  groups  of  birds  in  the  Great  Basin  to 
differences  in  vegetation  structure  and  composition  will  provide  information  on  how  species  that 
primarily  are  associated  with  shrublands,  coniferous  woodlands,  and  adjacent  riparian  areas  may 
respond  to  planned  removal  of  native  pinyon  ( Pinus  monophylla )  and  juniper  (primarily 
Juniperus  osteosperma  and  J.  occidental™ )  trees  on  as  many  as  120,000  km2  (30  million  acres). 

Future  research 

We  think  there  is  considerable  potential  to  conduct  novel,  empirical  tests  of  hypotheses  about  the 
information  content  of  occupancy  estimates.  As  noted  above,  theory  suggests  that  occupancy  is  a 
reliable  surrogate  measure  of  a  species’  abundance  and  that  abundance,  in  turn,  is  strongly 
related  to  probability  of  persistence.  However,  empirical  tests  of  the  hypothesized  relations 
between  measures  of  occurrence  and  persistence  have  not  been  conducted.  We  suggest  that  it 
would  be  valuable  and  feasible  to  test  whether  occupancy,  density  (abundance  per  unit  area),  and 
reproductive  success  are  correlated.  Given  that  maximizing  the  probability  of  species  persistence 
commonly  is  a  high  priority  for  resource  managers,  this  work  would  be  greatly  relevant  to 
decision-making.  Additionally,  theory  suggests  that  environmental  covariates  that  explain 
considerable  variation  in  observed  occupancy  may  be  interpreted  as  measures  of  habitat  quality, 
and  used  to  project  occupancy  in  unsampled  locations  or  time  periods.  Testing  this  hypothesis 
would  be  relevant  to  management  because  measurement  of  environmental  variables  typically  is 
more  feasible  than  measuring  occupancy  or  demography,  and  measures  of  environmental 
variables  often  have  less  uncertainty  than  measures  of  the  distribution,  survival,  or  reproduction 
of  individual  animals.  In  the  Chesapeake  Bay  Lowlands,  results  would  be  highly  relevant  to 
management  of  populations  of  white-tailed  deer,  which  have  substantial  effects  on  plant  and  bird 
communities.  In  any  part  of  the  Great  Basin,  our  results  would  be  applicable  to  planned 
management  treatments  on  public  lands,  such  as  removal  of  conifers  that  is  intended  to  increase 
the  probability  of  persistence  of  Greater  Sage-Grouse  ( Centrocercus  urophasianus )  and  benefit 
more  than  350  other  species  of  animals  (BLM  2015,  USFS  2015). 

We  initially  used  data  on  birds  and  butterflies  from  the  central  and  western  Great  Basin  to 
estimate  species  richness  of  multiple  taxonomic  groups  on  the  basis  of  a  single  group.  We  used 
these  data  because  they  best  supported  external  evaluation  of  the  spatial  and  temporal 
transferability  of  indicator  species.  Given  that  our  methods  are  fully  transferable,  although  the 
information  content  of  putative  indicator  species  likely  varies,  we  anticipate  applying  the 
methods  to  our  data  on  birds  and  butterflies  in  the  Chesapeake  Bay  Lowlands.  Moreover,  we 
anticipate  applying  multiple- state  occupancy  models  to  selected  species  of  birds  and  butterflies 
in  the  central  and  western  Great  Basin.  We  are  eager  to  explore  environmental  associations  with 
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abundances  of  these  taxonomic  groups  at  multiple  spatial  and  temporal  extents.  We  also  hope  to 
examine  whether  inferences  about  extents  associated  with  species  richness  and  beta  diversity  are 
consistent  with  inferences  about  extents  associated  with  abundance. 

We  also  think  it  would  be  feasible  and  worthwhile  to  explore  how  breeding  phenology, 
reproductive  success,  density,  and  occupancy  of  birds  in  the  Great  Basin  are  associated  with 
contemporary  and  projected  variation  in  weather  and  climate;  and  to  examine  how  the  interaction 
of  phenological  responses  to  climate  and  other  biological  responses  to  management-driven 
changes  in  land  cover  may  affect  the  feasibility  of  conserving  species  of  management  concern. 
For  example,  we  could  test  whether  our  data  on  subregional  trends  in  occupancy  and  abundance 
are  consistent  with  trends  across  the  Great  Basin  Bird  Conservation  Region,  and  whether  trends 
at  either  spatial  extent  are  associated  with  the  species’  biological  traits.  It  would  be 
straightforward  to  test  whether  subregional  or  Bird  Conservation  Region-level  trends  in 
occupancy  and  abundance  are  associated  with  published  climate-change  vulnerability  rankings. 
These  tests  would  be  of  interest  and  relevance  to  other  researchers  and  managers  given  recent 
development  of  methods  to  characterize  phenology  on  the  basis  of  occupancy  data.  They  would 
address  whether  local  responses  to  environmental  change  may  be  transferable  across  the  Great 
Basin  and  whether  responses  to  climate  change  potentially  can  be  predicted  on  the  basis  of  life 
history  information.  Furthermore,  we  believe  it  would  be  useful  to  test  whether  contemporary 
elevational  ranges  and  temporal  shifts  in  species’  elevational  ranges  are  associated  with  local 
climate,  and  whether  breeding  phenology  and  reproductive  success  are  associated  with  extremes, 
variability,  or  means  of  local  weather.  Moreover,  one  could  test  whether  the  nexus  between 
climate  change  and  proposed  vegetation  management  may  affect  the  breeding  phenology  or 
viability  of  species  of  concern. 

Although  the  Chesapeake  Bay  Lowlands  and  Great  Basin  are  ecologically  distinct,  some  species 
of  songbirds  in  either  region  have  similar  breeding  habitat  and  can  be  matched  phylogenetically. 
As  a  result,  one  can  examine  how  different  configurations  of  vegetation,  and  different  drivers  of 
changes  in  vegetation  structure  and  composition,  are  associated  with  population  dynamics  of 
similar  species.  For  example.  Hooded  Warbler  (Setophaga  citrina )  in  the  Chesapeake  Bay 
Lowlands  and  MacGillivray’s  Warbler  {Geothlypis  tolmiei )  in  the  Great  Basin  both  build  open- 
cup  nests  in  shrubs  and  are  associated  with  deciduous  woodlands.  As  another  illustration,  Wood 
Thrush  ( Hylocichla  mustelina )  in  the  Chesapeake  Bay  Lowlands  and  Cassin’s  Finch 
( Haemorhous  cassinii)  in  the  Great  Basin  both  build  cup  nests  fairly  high  in  trees  and  are 
associated  with  diverse  woodlands. 

Implementation 

We  regularly  shared  data,  inferences,  reports,  and  publications  with  resource  managers  and 
wildlife  biologists  at  Fort  Eustis,  Hawthorne  Army  Depot,  the  Marine  Corps  Mountain  Warfare 
Training  Center,  and  the  Naval  Facilities  Engineering  Command’s  Southwest  Central  Integrated 
Project  Team;  with  DoD’s  federal  partners,  including  the  US  Fish  and  Wildlife  Service,  Bureau 
of  Land  Management,  US  Forest  Service,  and  Great  Basin  Landscape  Conservation  Cooperative; 
and  with  state  and  local  agencies  and  citizens’  groups.  DoD  resource  managers  and  their 
commanding  officers  indicated  that  our  results  would  contribute  to  watershed  management 
programs  and  Integrated  Natural  Resources  Management  Plans.  We  also  responded  to  the 
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installations’  periodic  requests  for  information  on  topics  such  as  the  ecology  and  locations  where 
we  recorded  particular  species  of  rare  birds.  Additionally,  Fort  Eustis  was  interested  in  related 
work,  supported  by  other  sources,  on  relations  between  densities  of  ticks  and  white-tailed  deer. 

Our  results  may  inform  not  only  assessment  and  monitoring  of  native  species  but  land  use  and 
management  actions  intended  to  maintain  these  species.  For  example,  our  work  suggests  that  if 
maintenance  of  anurans  is  a  management  objective,  then  conservation  of  contemporary 
populations  in  situ  may  be  more  effective  than  restoration  or  construction  of  new  wetlands.  As 
another  illustration,  our  results  suggest  that  in  coastal  Virginia,  culling  of  white-tailed  deer  may 
increase  habitat  quality  for  species  of  birds  that  nest  in  shrubs  and  glean  insect  prey  from  foliage, 
including  several  species  of  conservation  concern.  We  found  that  it  may  be  useful  to  include  not 
only  larval  hostplants  but  plants  that  provide  nectar  for  adult  butterflies  in  assessment  and 
monitoring  programs  for  butterflies,  and  to  project  how  land  use  (e.g.,  road  maintenance,  seeding 
after  disturbances)  may  affect  these  aspects  of  butterfly  habitat.  More  generally,  environmental 
associations  with  species  richness  or  occupancy  of  anurans,  birds,  and  butterflies  that  were 
identified  in  our  work  can  be  used  to  project  how  natural  or  anthropogenic  phenomena  that  affect 
these  environmental  attributes  also  may  affect  native  species. 
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Appendix  A.  Supporting  data 

All  data  collected  in  the  course  of  this  project  are  being  archived. 

Data  collected  in  the  Chesapeake  Bay  Lowlands  will  be  deposited  in  the  College  of  William  and 
Mary’s  digital  archive  (https://digitalarchive.wm.eduj.  This  archive  is  the  College  of  William 
and  Mary  and  Virginia  Institute  of  Marine  Science’s  online  repository  of  research.  The  archive  is 
an  effort  to  collect,  preserve,  and  distribute  digital  material  related  to  and  produced  by  the 
university  and  its  students,  faculty,  and  staff,  and  is  accessible  to  the  public. 

Most  of  the  data  collected  in  the  Great  Basin  have  been  deposited  with  the  Forest  Service 
Research  Data  Archive  (see  Appendix  B).  Data  that  have  not  yet  been  deposited  will  be  added  to 
the  archive  in  the  form  of  new  editions  of  the  existing  data  packages. 
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